Tridiagonal matrices

A tridiagonal matrix is a matrix that has nonzero entries only on the main diagonal and on the adjacent off-diagonals. This special structure comes up frequently in applications. For example, the finite difference numerical solution to the heat equation leads to a tridiagonal system. Another application, the one we’ll look at in detail here, is natural cubic splines. And we’ll mention an interesting result with Fibonacci numbers in passing.

Because these matrices have a special structure, the corresponding linear systems are quick and easy to solve. Also, we can find a simple recurrence relation for their determinants.

Determinants

Let’s label the component of a tridiagonal matrix as below

where every component not shown is implicitly zero. We can expand the determinant of the matrix above using minors along the last row. This gives a recursive expression for the determinant

with initial conditions d₀ = 1 and d₋₁ = 0.

Note that if all the a‘s and b‘s are 1 and all the c‘s are −1, then you get the recurrence relation that defines the Fibonacci numbers. That is, the Fibonacci numbers are given by the determinant

Natural cubic splines

A cubic spline interpolates a set of data points with piecewise cubic polynomials. There’s a (potentially) different cubic polynomial over each interval between input values, all fitted together so that the resulting function, its derivative, and its second derivative are all continuous.

Suppose you have input points, called knots in this context, t₀, t₁, … t_n and output values y₀, y₁, … y_n. For the spline to interpolate the data, its value at t_i must be y_i.

A cubic spline then is a set of n cubic polynomials, one for each interval [t_i, t_i+1]. A cubic polynomial has four coefficients, so we have 4n coefficients in total. At each interior knot, t₁ through t_n−1, we have four constraints. Both polynomials that meet at t_i must take on the value y_i at that point, and the two polynomials must have the same first and second derivative at that point. That gives us 4(n − 1) equations. The value of the first polynomial is specified on the left end at t₀ the value of the last polynomial is specified at the right end at t_n. This gives us 4n − 2 equations in total.

We need two more equations. A clamped cubic spline specifies the derivatives at each end point. The natural cubic spline specifies instead that the second derivatives at each end are zero. What is natural about a natural cubic spline? In a certain sense it is the smoothest curve interpolating the specified points. With these boundary conditions we now have as many constraints as degrees of freedom.

So how would we go about finding the coefficients of each polynomial? Our task will be much easier if we parameterize the polynomials cleverly to start with. Instead of powers of x, we want powers of (x − t_i) and (t_i+1 − x) because these expressions are 0 on different ends of the interval [t_i, t_i+1]. It turns out we parameterize the spline over the ith interval as

where h_i = t_i+1 − t_i, the length of the ith interval.

This may seem unmotivated, and no doubt it is cleaner than the first thing someone probably tried, but it’s the kind of thing you’re lead to when you try to make the derivation work out smoothly.

The basic form is  powers of (x − t_i) and (t_i+1 − x), each to the first and third powers, for reasons given above. Why the 6’s in the denominators? They’re not strictly necessary, but they cancel out when we take second derivatives. Let’s look at the second derivative.

Note how when we stick in t_i the first term is zero and the second is z_i, and when we stick in t_i+1 the first term is z_i+1 and the second is zero.

We can now write down the system of equations for the z‘s. We have z₀ = z_n from the natural cubic spline condition, and for 1 ≤ i ≤ n − 1 we have

where

Note that this is a tridiagonal system because the ith equation only involves z‘s with subscripts i − 1, i, and i + 1.

Because of its tridiagonal structure, these equations can be solved simply and efficiently, much more efficiently than a general system of equations.

Related math posts

• Applied linear algebra
• Chebyshev interpolation
• Help with interpolation

Fitting a polynomial to a function at more points might not produce a better approximation. This is Faber’s theorem, something I wrote about the other day.

If the function you’re interpolating is smooth, then interpolating at more points may or may not improve the fit of the interpolation, depending on where you put the points. The famous example of Runge shows that interpolating

f(x) = 1 / (1 + x²)

at more points can make the fit worse. When interpolating at 16 evenly spaced points, the behavior is wild at the ends of the interval.

Here’s the Python code that produced the plot.

    import matplotlib.pyplot as plt
    from scipy import interpolate, linspace

    def cauchy(x):
        return (1 + x**2)**-1

    n = 16
    x = linspace(-5, 5, n) # points to interpolate at

    y = cauchy(x)
    f = interpolate.BarycentricInterpolator(x, y)

    xnew = linspace(-5, 5, 200) # points to plot at
    ynew = f(xnew)
    plt.plot(x, y, 'o', xnew, ynew, '-')
    plt.show()

However, for smooth functions interpolating at more points does improve the fit if you interpolate at the roots of Chebyshev polynomials. As you interpolate at the roots of higher and higher degree Chebyshev polynomials, the interpolants converge to the function being interpolated. The plot below shows how interpolating at the roots of T₁₆, the 16th Chebyshev polynomial, eliminates the bad behavior at the ends.

To make this plot, we replaced x above with the roots of T₁₆, rescaled from the interval [-1, 1] to the interval [-5, 5] to match the example above.

    x = [cos(pi*(2*k-1)/(2*n)) for k in range(1, n+1)]
    x = 5*array(x)

What if the function we’re interpolating isn’t smooth? If the function has a step discontinuity, we can see Gibbs phenomena, similar to what we saw in the previous post. Here’s the result of interpolating the indicator function of the interval [−1, 1] at 100 Chebyshev points. We get the same “bat ears” as before.

Related: Help with interpolation

I just got an evaluation copy of The Best Writing on Mathematics 2017. My favorite chapter was Inverse Yogiisms by Lloyd N. Trefethen.

Trefethen gives several famous Yogi Berra quotes and concludes that

Yogiisms are statements that, if taken literally, are meaningless or contradictory or nonsensical or tautological—yet nevertheless convey something true.

An inverse yogiism is the opposite,

[a] statement that is literally true, yet conveys something false.

What a great way to frame a chapter! Now that I’ve heard the phrase, I’m trying to think of inverse yogiisms. Nothing particular has come to mind yet, but I feel like there must be lots of things that fit that description. Trefethen comes up with three inverse yogiisms, and my favorite is the middle one: Faber’s theorem on polynomial interpolation.

Faber’s theorem is a non-convergence result for interpolants of continuous functions. Trefethen quotes several numerical analysis textbooks that comment on Faber’s theorem in a way that implies an overly pessimistic interpretation. Faber’s theorem is true for continuous functions in general, but if the function f  being interpolated is smooth, or even just Lipschitz continuous, the theorem doesn’t hold. In particular, Chebyshev interpolation produces a sequence of polynomials converging to f.

A few years ago I wrote a blog post that shows a famous example due to Carle Runge that if you interpolate f(x) = 1/(1 + x²) over [−5, 5] with evenly spaced nodes, the sequence of interpolating polynomials diverges. In other words, adding more interpolation points makes the fit worse.

Here’s the result of fitting a 16th degree polynomial to f  at evenly spaced nodes.

The error near the ends is terrible, though the fit does improve in the middle. If instead of using evenly spaced nodes you use the roots of Chebyshev polynomials, the interpolating polynomials do in fact converge, and converge quickly. If the kth derivative of f has bounded variation, then the error in interpolating f at n points is O(n^−k).

Say you have a function f(x) and you want to find a polynomial p(x) that agrees with f(x) at several points. Given a set of points

x₀, x₁, x₂, …, x_n

you can always find a polynomial of degree n such that p(x_i) = f(x_i) for i = 0, 1, 2, …, n. It seems reasonable that the more points you pick, the better the interpolating polynomial p(x) will match the function f(x). If the two functions match exactly at a lot of points, they should match well everywhere. Sometimes that’s true and sometimes it’s not.

Here is a famous example due to Carle Runge. Let f(x) = 1/(1 + x²) and let p_n be the polynomial that interpolates f(x) at n+1 evenly spaced nodes in the interval [−5, 5]. As n becomes larger, the fit becomes worse.

Here’s a graph of f(x) and p₉(x). The smooth blue line is f(x) and the wiggly red line is p₉(x).

Here’s the analogous graph for p₁₆(x).

The fit is improving in the middle. In fact, the curves agree to within the thickness of the plot line from say −1 to 1. But the fit is so bad in the tails that the graph had to be cut off. Here’s another plot of f(x) and p₁₆(x) on a different scale to show how far negative the polynomial dips.

The problem is the spacing of the nodes. Interpolation errors are bad for evenly spaced nodes.

Update: This post explains in a little more depth why this particular function has problems and gives another example where interpolation at evenly-spaced nodes behaves badly.

If we interpolate f(x) at different points, at the Chebyshev nodes, then the fit is good.

The Chebyshev nodes on [−1, 1] are x_i = cos( π i / n ). Here we multiplied these nodes by 5 to scale to the interval [−-5, 5].

If the function f(x) is absolutely continuous, as in our example, then the interpolating polynomials converge uniformly when you interpolate at Chebyshev nodes. However, ordinary continuity is not enough. Given any sequence of nodes, there exists a continuous function such that the polynomial interpolation error grows like log(n) as the number of nodes n increases.

Some numerical integration methods are based on interpolating polynomials: fit a polynomial to the integrand, then integrate the polynomial exactly to approximate the original integral. The examples above suggest that increasing the order of such integration methods might not improve accuracy and might even make things worse.

Related: Need help with interpolation?