Polynomial interpolation in Haskell

A classical problem in scientific computing is the following: Given data points \((x_0,y_0),\ldots,(x_n,y_n)\), obtain a polynomial \(p\) such that \(p(x_i) = y_i\) for \(i = 0,\ldots,n\). We say that \(p\) interpolates the data points and that \(p\) is an interpolating polynomial or simply an interpolant.

Often, the \(y_i\)’s are obtained by evaluating a function that is difficult to compute or describe succintly at sample input values (the \(x_i\)’s in this case). It is hoped that a polynomial that agree in values with the unknown function at the chosen inputs gives a reasonably good approximation of the function so that further analyses involving the function can also be approximated.

We assume that the \(x_i\)’s are distinct. Since there are \(n+1\) points, we do not need to look at polynomials of degree more than \(n\). In particular, we seek values \(a_0,\ldots,a_n\) such that \(p(x) = a_0 + a_1x + \ldots + a_n x^n\) satisfies \(p(x_i) = y_i\) for \(i = 0,\ldots,n\), which can be written in matrix form as \[\begin{equation*} \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix} \end{equation*}\]

The coefficient matrix is well known as a Vandermonde matrix with determinant \(\displaystyle\prod_{0 \leq i < j \leq n} (x_j-x_i)\) which is nonzero when the \(x_i\)’s are distinct. Hence, the coefficient matrix is invertible and the matrix equation has a unique solution. It follows that there is a unique polynomial of degree at most \(n\) that interpolates the given data points. For example, say we have the following data points:

\(i\)	\((x_i,y_i)\)
\(0\)	\((2,1)\)
\(1\)	\((5,10)\)
\(2\)	\((7,-24)\)
\(3\)	\((8,-17)\)

We have the following system of equations: \[\begin{equation*} \begin{bmatrix} 1 & 2 & 2^2 & 2^3 \\ 1 & 5 & 5^2 & 5^3 \\ 1 & 7 & 7^2 & 7^3 \\ 1 & 8 & 8^2 & 8^3 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_3 \\ a_4 \end{bmatrix} = \begin{bmatrix} 1 \\ 10 \\ -24 \\ -17 \end{bmatrix}. \end{equation*}\] Solving gives \(a_0 = -185\), \(a_1 = 149\), \(a_2 = -32\), \(a_3 = 2\). Hence, we have the interpolating polynomial \[-185 + 149x - 32x^2 + 2x^3.\]

In general, solving the matrix equation directly (as a system of linear equations) is not necessarily a good idea because the matrix is notorious for being ill-conditioned. In this article, we will look at polynomial interpolation using the so-called Newton form which has the advantage of being able to interpolate in an online manner; that is, construct the polynomial as the data points arrive one by one. We will then implement the contruction in Haskell.

Choosing a good basis

Let \(\mathcal{P}_n\) denote the set of polynomials in \(x\) of degree at most \(n\) with real coefficients. It is known that \(\mathcal{P}_n\) form a vector space over the real numbers \(\mathbb{R}\) though we do not need this fact here. What we need is that for any given set of polynomials \(p_0(x),\ldots, p_n(x)\) where \(p_i\) is of degree \(i\), we can find \(c_0,\ldots,c_n \in \mathbb{R}\) so that \[p = c_0p_0 + c_2p_2 + \cdots + c_np_n\] for any given \(p \in \mathcal{P}_n\). Observe that it suffices to focus on \(p = 1,x, \ldots, x^n\) in order to prove this statement.

The basis case is when \(p = 1\). We can set \(c_0 = 1/p_0\) and \(c_i = 0\) for \(i = 1,\ldots, n\).

Consider \(p = x^k\) where \(k \geq 1\). Write \(p_k(x)\) as \(b_kx^k + q(x)\) where \(q\) has degree at most \(k-1\). By the induction hypothesis, there exist \(c'_0,\ldots,c'_{k-1} \in \mathbb{R}\) so that \[q = c'_0p_0 + \cdots + c'_{k-1}p_{k-1}.\] Hence, \[\begin{align*} x^k & = \frac{1}{b_k}\left(p_k(x) - q(x) \right) \\ & = \frac{1}{b_k}p_k(x) - \frac{1}{c'_{k-1}}p_{k-1}(x) - \cdots - \frac{c'_0}{b_0}p_0(x) \\ & = c_k x^k + c_{k-1} x^{k-1} + \cdots + c_1 x + c_0 \end{align*}\] where \(c_i = -\frac{c'_i}{b_k}\) for \(i = 0,\ldots,k-1\) and \(c_k = \frac{1}{b_k}\). This completes the induction.

Using this result, we see that the interpolating polynomial of the data points \((x_0,y_0),\ldots,(x_n,y_n)\) can be written in the so-called Newton form \[c_0p_0(x) + c_1p_1(x) + c_2p_2(x) + \cdots + c_np_n(x)\] where \(p_0(x) = 1\) and \(p_i(x) = (x-x_0)(x-x_1)\cdots(x-x_{i-1})\) for \(i = 1,\ldots,n\).

An attractive feature of this form is that \(p_i(x_j) = 0\) for \(j \geq i\), implying that \[c_0p_0(x) + c_1p_1(x) + \cdots + c_kp_k(x)\] is the interpolating polynomial for the data points \((x_0,y_0),\ldots,(x_k,y_k)\)!

In other words, we can compute the scalars \(c_0,\ldots,c_n\) one at a time as each data point comes in. In particular, if \(c_0,\ldots,c_k\) have already been computed from the data points \((x_0,y_0),\ldots,(x_k,y_k)\), we can obtain \(c_{k+1}\) by computing

Example

With the data points in the example above, we have \[\begin{align*} c_0 & = y_0 = 1, \\ c_1 & = \frac{y_1 - c_0}{x_1-x_0} = \frac{10 - 1}{5-2} = 3, \\ c_2 & = \frac{y_2 - (c_0+c_1(x_2-x_0))}{(x_2-x_0)(x_2-x_1)} = \frac{-24 - 16}{10} = -4, \\ c_3 & = \frac{y_3-(c_0+(x_3-x_0)(c_1+c_2(x_3-x_1)))}{(x_3-x_0)(x_3-x_1)(x_3-x_2)} = \frac{-17 + 53}{18} = 2. \end{align*}\]

Hence, the polynomial \[\begin{align*} p(x) & = 1 + (x - 2)(3 + (x-5)(-4 + 2(x-7))) \\ & = 1 + 3(x-2) - 4(x-2)(x-5) + 2(x-2)(x-5)(x-7) \end{align*}\] interpolates the given data points.

Implementation in Haskell

{-# LANGUAGE StrictData #-}

module PolynomialInterpolation where

import Data.Vector (Vector)
import qualified Data.Vector as V
import Data.List (foldl')

data NewtonForm a =   Empty 
                    | NewtonForm { nf_cs :: !(Vector a), nf_xs :: !(Vector a) }
                    deriving (Eq, Show)

We first develop a function that evaluates the Newton form at a given input. The input will be a list consisting of the pairs \((c_0,x_0), \ldots, (c_n,x_n)\). (Note that \(x_n\) will not be used in the evaluation.) We use Vector for holding such a list.

evalNewton :: Fractional a => NewtonForm a -> a -> a
evalNewton Empty _ = 0
evalNewton nf x =
  case V.length cs of
    1 -> c
    _ -> c + (x - x') * evalNewton (NewtonForm (V.tail cs) (V.tail xs)) x
  where
    cs = nf_cs nf
    xs = nf_xs nf
    c = V.head cs
    x' = V.head xs

For the recursive step, we develop a function that takes a current interpolant in Newton form a new data point, and outputs the new Newton form:

nextNewtonForm :: Fractional a =>    NewtonForm a -- the current interpolant
                                  -> (a, a) -- next data point
                                  -> NewtonForm a
nextNewtonForm Empty (x,y) = NewtonForm (V.singleton y) (V.singleton x)
nextNewtonForm nf (x,y) = 
    NewtonForm (V.snoc cs c) (V.snoc xs x)
  where
    cs = nf_cs nf
    xs = nf_xs nf
    denominator = V.product (V.map (x - ) xs)
    numerator = y - evalNewton nf x
    c = numerator / denominator

A simple fold will then give us a function that computes the interpolant in Newton form from a list of data points:

newtonInterpolation :: Fractional a => [(a, a)] -> NewtonForm a
newtonInterpolation = foldl' nextNewtonForm Empty

λ> :r
λ> newtonInterpolation [(2,1),(5,10),(7,-24),(8,-17)]
NewtonForm {nf_cs = [1.0,3.0,-4.0,2.0], nf_xs = [2.0,5.0,7.0,8.0]}

The coefficients are \(1,3,-4,2\), which agree with what we had obtained before.