Using Chebyshev points for Lagrange interpolation

Lagrange interpolation is a useful polynomial interpolation formula that lets us fit to an arbitrary function. While we might expect that increasing the number of points would decrease the average error of our polynomial approximation, this intuition fails because of Runge's phenomenon. If we pick our points carefully, though, we can overcome this phenomenon and still achieve a great approximation.

Specifically, in this post, we're going to approximate the function \(\frac{1}{1 + x^6}\) over \([0, 5]\) by using two different sets of points: a set of equidistant points across the interval vs. Chebyshev nodes. We'll compare the average error (MSE) and the maximum error over the interval.

Helper functions

We'll add a function that returns the Lagrange interpolating polynomial for a set of points:

from operator import mul
from functools import reduce

def lagrange(pts, f):
    """ Return a Lagrange interpolating polynomial for f
    at the points pts. """

    def basis(others):
        # w_{i}(x) = (x - x_{1}) * ... * (x - x_{i-1}) * (x - x_{i+1}) * ...
        return lambda x : reduce(mul, [x - o for o in others])

    # Get each basis polynomial and coefficients
    bases = [basis(pts[:i] + pts[i + 1:]) for i in range(len(pts))]
    # a_{i} = f(x_{i}) / w_{i}(x_{i})
    weights = [f(pt) / bases[i](pt) for i, pt in enumerate(pts)]

    return lambda x : sum([weights[i] * bases[i](x) for i in range(len(pts))])

This is a basic implementation of Lagrange interpolation - nothing special here.

Next, we'll add a couple of helper functions to compute the average and maximum error:

def average_error(p, f):
    # Get the MSE of p over 1000 points on the interval [0, 5]
    N = 1000
    error = 0
    for pt in [5 * i / N for i in range(N)]:
        error += (p(pt) - f(pt))**2
    return error / N

def max_error(p, f):
    error = float("-inf")
    for pt in [5 * i / 1000 for i in range(1000)]:
        curr_error = abs(p(pt) - f(pt))
        if curr_error > error:
            error = curr_error
    return curr_error

Computing for equidistant vs. Chebyshev points

Finally, we'll add our function to approximate:

def g(x):
    return 1 / (1 + x**6)

And some functions to compute and evaluate the interpolating polynomial for N points, whether equidistant or Chebyshev:

def evaluate(pts, f, tag):
    print("Interpolating for {0}".format(tag))
    poly = lagrange(pts, f)
    print("Average error: {0}".format(average_error(poly, f)))
    print("Max error: {0}".format(max_error(poly, f)))

def equidistant(n):
    # Interpolate for n points from [0, 5]
    pts = [5 * i / n for i in range(n)]
    evaluate(pts, g, "{0} equidistant points".format(n))

import math
def chebyshev(n):
    # Interpolate for n points from [0, 5]
    pts = [(5 + 5 * math.cos((2 * i + 1) * math.pi / (2 * n + 2))) / 2 for i in range(n)]
    evaluate(pts, g, "{0} Chebyshev points".format(n))

n_values = [3, 5, 10, 20, 30]
for n in n_values:
    equidistant(n)
    print("")
    chebyshev(n)
    print("")

Evaluation

Here's what the results show:

Mean Squared Error:

Number of points	Equidistant	Chebyshev
3	0.083449016830	0.123504372195
5	0.082739268475	0.08960932097
10	0.64207073475	0.0127952042780
20	26.4328328320	0.0005982536404
30	435.529425177	9.3148785757e-06

Maximum Error:

Number of points	Equidistant	Chebyshev
3	0.86442582784	0.0030308889852
5	1.34512664557	0.0063971895249
10	6.2268209092	0.000121486335635
20	61.205287745	2.1722048813e-05
30	312.690373601	4.1559874016e-06

I initially wanted to graph these results, but as you can see, the error for the Chebyshev nodes is so low that it would hardly be visible.

All in all, clearly it's preferable to use Chebyshev nodes. This is especially clear as the number of points increases. I hope this was helpful!