Cubic spline interpolation exercise

Large language models (predictably) learn to represent the semantic meaning of sentences.

Cubic spline interpolation exercise

I somehow stumbled upon cubic spline interpolation quite accidentally, but it caught my interest with all of the neat linear algebra and scientific computing techniques behind it. This got me excited to revise materials from my undergrad classes and try to implement the method in Python from scratch.

Resources:

1. The scipy scipy implementation of the method.
2. Cubic Spline Interpolation on Wikiversity.

Imports

Toggle cells below if you want to see what imports are being made.

%load_ext autoreload
%autoreload 2

import numpy as np
import plotly.graph_objects as go
from scipy.interpolate import CubicSpline

# Ensures we can render plotly plots with quarto
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook_connected"

Data

Cubic spline interpolation with scipy

x_train = np.linspace(0, 10, 11)
y_train = np.sin(x_train) + np.random.normal(0, 0.5, 11)

cs = CubicSpline(x_train, y_train, bc_type="natural", extrapolate=True)

x_test = np.linspace(-1, 11, 121)
y_test = cs(x_test)

# Create a figure
fig = go.Figure()

# Add scatter plot
fig.add_trace(go.Scatter(x=x_train, y=y_train, mode='markers', name='Original points'))

# Add line plot
fig.add_trace(go.Scatter(x=x_test, y=y_test, mode='lines', name='Fitted line'))

# Add titles and labels
fig.update_layout(title="Scipy interpolation", xaxis_title="x", yaxis_title="y")

# Show the plot
fig.show()

Custom implementation of cubic spline interpolation

def tridiagonal_linear_system_solver(d_lower, d_main, d_upper, b):
    """
    Solve a tridiagonal linear system given the main, lower, and upper diagonals, as well as the vector b. The Thomas
    algorithm is used.
    """

    n = len(d_main)

    # Forward sweep
    for i in range(1, n):
        w = d_lower[i - 1] / d_main[i - 1]
        d_main[i] -= w * d_upper[i - 1]
        b[i] -= w * b[i - 1]

    # Back substitution
    x = np.zeros(n)
    x[-1] = b[-1] / d_main[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (b[i] - d_upper[i] * x[i + 1]) / d_main[i]

    return x


class CubicSplineCustom:

    def __init__(self, x, y):
        if not np.all(np.diff(x) > 0):
            raise ValueError("x values must be in ascending order.")

        self.x = np.array(x)
        self.y = np.array(y)
        self._find_coeffs()

    def _find_coeffs(self):
        # Find the relevant diagonals and from self.x and self.y assuming natural boundary conditions
        d_lower, d_main, d_upper, b, h = self._compute_diagonals_and_b()

        # Compute the M_i's for i = 1, n-1 (since M_0 and M_n are assumed to be 0). This will require solving a
        # tridiagonal linear system
        ms = tridiagonal_linear_system_solver(d_lower, d_main, d_upper, b)

        # Compute the coefficients of the cubic polynomials
        self.c = self._compute_coefficients_from_second_derivatives(ms, h)

    def _compute_diagonals_and_b(self):
        x, y = self.x, self.y
        h = np.diff(x)

        # Compute the diagonals for the tridiagonal matrix
        d_lower = h.copy()
        d_lower[-1] = 0  # one of the naturals BCs
        d_upper = h.copy()
        d_upper[0] = 0

        d_main = 2 * np.ones(len(x))
        d_main[1:-1] *= h[:-1] + h[1:]

        b = np.zeros(len(x))
        y_diff = np.diff(y)
        b[1:-1] = 6 * (y_diff[1:] / h[1:] - y_diff[:-1] / h[:-1])

        return d_lower, d_main, d_upper, b, h

    def _compute_coefficients_from_second_derivatives(self, ms, h):
        coeffs = np.zeros((4, len(self.x) - 1))
        coeffs[0, :] = (ms[1:] - ms[:-1]) / (6 * h)
        coeffs[1, :] = ms[:-1] / 2
        coeffs[2, :] = (self.y[1:] - self.y[:-1]) / h - (ms[1:] + 2 * ms[:-1]) * h / 6
        coeffs[3, :] = self.y[:-1]

        return coeffs

    def _get_index(self, x):
        """Performs binary search"""
        low, high = 0, len(self.x) - 1

        while low < high:
            mid = (low + high) // 2
            if self.x[mid] <= x:
                low = mid + 1
            else:
                high = mid

        return min(max(low - 1, 0), len(self.x) - 2)

    def interpolate(self, x: float):
        # Find which polynomial is appropriate and evaluate at x
        idx = self._get_index(x)
        dx = x - self.x[idx]
        return (self.c[0, idx] * dx + self.c[1, idx]) * dx**2 + self.c[2, idx] * dx + self.c[3, idx]

custom_spline = CubicSplineCustom(x_train, y_train)

assert cs.c.shape == custom_spline.c.shape
assert np.allclose(cs.c, custom_spline.c)
np.sqrt(((cs.c - custom_spline.c) ** 2).mean())

2.4075003542614415e-16

assert np.allclose(custom_spline.interpolate(-1), cs(-1))
assert np.allclose(custom_spline.interpolate(86500), cs(86500))

# Create a figure
fig = go.Figure()

# Add scatter plot
fig.add_trace(go.Scatter(x=x_train, y=y_train, mode="markers", name="Original data"))

# Add line plot
fig.add_trace(go.Scatter(x=x_test, y=y_test, mode="lines", name="Fitted line"))

y_test_custom = [custom_spline.interpolate(x) for x in x_test]
assert np.allclose(y_test, y_test_custom)
fig.add_trace(
    go.Scatter(x=x_test, y=y_test_custom, mode="lines", name="Fitted line (custom)", line=dict(dash="longdash"))
)

# Add titles and labels
fig.update_layout(title="Scipy and custom interpolations", xaxis_title="x", yaxis_title="y")

# Show the plot
fig.show()