Spline Regression › Formula

The spline regression model extends linear regression by using piecewise polynomial basis functions instead of simple polynomial terms. We'll start with the most intuitive form and progressively build to the complete mathematical framework, explaining every piece of notation along the way.

Let's begin by understanding what a spline function is. A spline function $S(x)$ is a special type of function that is constructed by joining together multiple polynomial pieces. Think of it like connecting different curved segments to create one smooth, continuous curve.

The mathematical definition of a spline function $S(x)$ of degree $d$ with knots at positions $t_1, t_2, \ldots, t_k$ is:

$$S(x) = \begin{cases} P_1(x) & \text{if } x \leq t_1 \\ P_2(x) & \text{if } t_1 < x \leq t_2 \\ \vdots & \vdots \\ P_{k+1}(x) & \text{if } x > t_k \end{cases}$$

Let's break down every piece of this notation:

• $S(x)$: This is our spline function. It takes an input value $x$ and returns a corresponding output value. The function is defined piecewise, meaning it has different formulas for different ranges of $x$.

• $d$: This represents the degree of the spline. If $d = 1$, we have linear splines (straight line segments). If $d = 2$, we have quadratic splines (parabolic segments). If $d = 3$, we have cubic splines (cubic polynomial segments). Cubic splines are most commonly used because they provide a good balance between flexibility and smoothness.

• $t_1, t_2, \ldots, t_k$: These are called knots or knot points. They are specific values of $x$ where we transition from one polynomial piece to another. Think of them as "break points" where the function changes its behavior. For example, if $t_1 = 3$ and $t_2 = 7$, then:

  • For all $x$ values less than or equal to 3, we use polynomial $P_1(x)$
  • For all $x$ values between 3 and 7, we use polynomial $P_2(x)$
  • For all $x$ values greater than 7, we use polynomial $P_3(x)$

• $P_i(x)$: Each $P_i(x)$ is a polynomial of degree $d$ or less. For example, if $d = 3$ (cubic splines), then each $P_i(x)$ might look like:

  $$P_i(x) = a_i + b_i x + c_i x^2 + d_i x^3$$

  where $a_i, b_i, c_i, d_i$ are coefficients that we need to determine.

Interval Notation:

The notation $x \leq t_1$, $t_1 < x \leq t_2$, and $x > t_k$ uses mathematical interval notation:

• $x \leq t_1$: This means "all values of $x$ that are less than or equal to $t_1$"
• $t_1 < x \leq t_2$: This means "all values of $x$ that are greater than $t_1$ AND less than or equal to $t_2$"
• $x > t_k$: This means "all values of $x$ that are greater than $t_k$"

The key insight is that we're not fitting a single polynomial across the entire range, but rather different polynomials in different regions, connected at the knot points. This allows us to capture local patterns in the data while maintaining smoothness across the entire function.

This visualization demonstrates the fundamental concept of spline regression: instead of using a single complex polynomial across the entire range, we use simpler polynomials in different regions. Each colored segment represents a different cubic polynomial, and these pieces connect smoothly at the knots (red dashed lines). The complete spline function is the combination of all these pieces, creating a flexible curve that can adapt to local patterns in the data.

The most common form is cubic spline regression, where we use cubic polynomials in each segment. The cubic spline regression model can be written as:

$$y = \beta_0 + \sum_{j=1}^{p} \beta_j B_j(x) + \epsilon$$

where:

• $y$: the response variable (the value we want to predict)
• $\beta_0$: the intercept term (baseline value when all basis functions are zero)
• $\beta_j$: the coefficient for the $j$-th B-spline basis function (weight assigned to basis function $j$)
• $B_j(x)$: the $j$-th B-spline basis function evaluated at input $x$ (building block of the spline)
• $p$: the total number of B-spline basis functions (determined by the number of knots and degree)
• $x$: the input variable or predictor
• $\epsilon$: the error term (captures unexplained variation in the response)

The basis functions $B_j(x)$ are constructed to ensure continuity and smoothness at the knot points.

The most elegant way to construct spline basis functions is through B-splines (Basis splines). B-splines are special functions that serve as the "building blocks" for constructing any spline function. Think of them as standardized, reusable components that we can combine in different ways to create the spline we want.

Why B-splines?

Instead of directly working with the piecewise polynomials $P_i(x)$, B-splines provide a more systematic approach. We can express any spline function as a linear combination of B-spline basis functions:

$$S(x) = \sum_{j=1}^{p} \beta_j B_j(x)$$

where:

• $S(x)$: the spline function evaluated at input value $x$ (the complete fitted curve)
• $B_j(x)$: the $j$-th B-spline basis function evaluated at $x$ (building block function)
• $\beta_j$: the coefficient for the $j$-th basis function (weight determined through regression)
• $p$: the total number of B-spline basis functions used in the model
• $j$: index ranging from 1 to $p$ (identifies which basis function)

Recursive Construction:

B-splines are constructed recursively, starting from simple step functions and building up to more complex, smooth functions. Let's understand this step by step.

Zero-degree B-splines (Step Functions):

The foundation of B-splines starts with zero-degree B-splines, which are essentially step functions:

$$B_{i,0}(x) = \begin{cases} 1 & \text{if } t_i \leq x < t_{i+1} \\ 0 & \text{otherwise} \end{cases}$$

Let's break down this notation:

• $B_{i,0}(x)$: This represents the $i$-th zero-degree B-spline function. The subscript $i$ tells us which knot interval this function is associated with, and the $0$ indicates it's a zero-degree (constant) function.

• $t_i \leq x < t_{i+1}$: This is the interval where the function equals 1. Notice the different inequality signs:

  • $\leq$ (less than or equal to) on the left means the interval includes the left endpoint $t_i$
  • $<$ (less than) on the right means the interval excludes the right endpoint $t_{i+1}$

• "otherwise": This means the function equals 0 everywhere else.

Example of Zero-degree B-splines:

If we have knots at $t_0 = 0, t_1 = 2, t_2 = 5, t_3 = 8$, then:

• $B_{0,0}(x) = 1$ when $0 \leq x < 2$, and $0$ elsewhere
• $B_{1,0}(x) = 1$ when $2 \leq x < 5$, and $0$ elsewhere
• $B_{2,0}(x) = 1$ when $5 \leq x < 8$, and $0$ elsewhere

This visualization demonstrates why cubic splines are the most commonly used: they provide C² continuity, meaning the function itself, its first derivative (slope), and its second derivative (curvature) are all continuous at the knots. Linear splines (C⁰) show visible corners at knots where the slope changes abruptly. Quadratic splines (C¹) have smooth slopes but visible changes in curvature. Cubic splines (C²) provide the smoothest transitions with no visible discontinuities in the function, slope, or curvature. This smoothness makes cubic splines ideal for modeling natural phenomena and physical relationships.

Now that we understand how B-splines work, let's see how we can express the spline regression model in matrix form. This will show us how to compute spline regression using linear algebra.

The spline regression model can be written in matrix form as:

$$\mathbf{y} = \mathbf{B}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$

where:

• $\mathbf{y}$: response vector of dimension $n \times 1$ (observed values we want to predict)
• $\mathbf{B}$: spline basis matrix of dimension $n \times p$ (evaluated basis functions at all data points)
• $\boldsymbol{\beta}$: coefficient vector of dimension $p \times 1$ (weights for each basis function)
• $\boldsymbol{\epsilon}$: error vector of dimension $n \times 1$ (unexplained variation in responses)
• $n$: number of observations in the dataset
• $p$: number of B-spline basis functions

This equation looks very similar to the standard linear regression equation, but with an important difference: instead of using simple polynomial terms, we're using B-spline basis functions. Let's break down each component:

The Response Vector:

A column vector containing all our observed response values:

$$\mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}$$

where $y_i$ is the observed response value for the $i$-th observation.

The Spline Basis Matrix:

$$\mathbf{B} = \begin{bmatrix} B_1(x_1) & B_2(x_1) & \cdots & B_p(x_1) \\ B_1(x_2) & B_2(x_2) & \cdots & B_p(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ B_1(x_n) & B_2(x_n) & \cdots & B_p(x_n) \end{bmatrix}$$

This is the most important part of the equation. The matrix $\mathbf{B}$ contains the values of all B-spline basis functions evaluated at all our data points:

• Rows: Each row corresponds to one observation. Row $i$ contains the values of all basis functions evaluated at $x_i$.
• Columns: Each column corresponds to one B-spline basis function. Column $j$ contains the values of $B_j(x)$ evaluated at all data points.
• Element $B_j(x_i)$: The entry in row $i$, column $j$ represents the value of the $j$-th B-spline basis function evaluated at the $i$-th data point.

Example of the Basis Matrix:

If we have 3 observations at $x_1 = 1, x_2 = 3, x_3 = 5$ and 4 B-spline basis functions, our basis matrix might look like:

$$\mathbf{B} = \begin{bmatrix} B_1(1) & B_2(1) & B_3(1) & B_4(1) \\ B_1(3) & B_2(3) & B_3(3) & B_4(3) \\ B_1(5) & B_2(5) & B_3(5) & B_4(5) \end{bmatrix} = \begin{bmatrix} 0.8 & 0.2 & 0.0 & 0.0 \\ 0.1 & 0.7 & 0.2 & 0.0 \\ 0.0 & 0.2 & 0.6 & 0.2 \end{bmatrix}$$

Notice how each row sums to approximately 1 (due to the partition of unity property), and most entries are 0 (due to the compact support property).

The Coefficient Vector:

$$\boldsymbol{\beta} = \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{bmatrix}$$

This vector contains the coefficients we need to determine through regression, where $\beta_j$ is the weight assigned to the $j$-th B-spline basis function.

The Error Vector:

$$\boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{bmatrix}$$

This vector contains the error terms (also called residuals), where $\epsilon_i$ represents the difference between the observed value $y_i$ and the predicted value for observation $i$.

Matrix Multiplication:

The expression $\mathbf{B}\boldsymbol{\beta}$ represents matrix multiplication. For each observation $i$, this gives us:

$$(\mathbf{B}\boldsymbol{\beta})_i = \sum_{j=1}^{p} B_j(x_i)\beta_j = B_1(x_i)\beta_1 + B_2(x_i)\beta_2 + \cdots + B_p(x_i)\beta_p$$

where:

• $(\mathbf{B}\boldsymbol{\beta})_i$: the predicted value for observation $i$ (the $i$-th element of the product)
• $B_j(x_i)$: the value of the $j$-th basis function at the $i$-th data point (element from row $i$, column $j$ of $\mathbf{B}$)
• $\beta_j$: the coefficient for the $j$-th basis function (element $j$ of $\boldsymbol{\beta}$)

This is exactly our spline function evaluated at $x_i$:

$$S(x_i) = \sum_{j=1}^{p} \beta_j B_j(x_i)$$

So the matrix equation $\mathbf{y} = \mathbf{B}\boldsymbol{\beta} + \boldsymbol{\epsilon}$ is saying that for each observation $i$:

$$y_i = S(x_i) + \epsilon_i$$

In other words, the observed value equals the spline function value plus some random error.

For better behavior at the boundaries (the edges of our data), we often use natural cubic splines. These impose additional constraints that help prevent the spline from behaving erratically when we try to predict values outside the range of our training data.

The Problem with Regular Cubic Splines: Regular cubic splines can sometimes extrapolate (predict beyond the data range) in ways that don't make physical sense. For example, if we're modeling temperature over time, a regular cubic spline might predict temperatures that are impossibly high or low when we try to forecast into the future.

The Natural Cubic Spline Solution: Natural cubic splines solve this problem by imposing boundary conditions that make the spline behave more reasonably at the edges. Specifically, natural cubic splines require that the second derivative is zero at the boundary knots.

Derivatives:

Before we can understand the constraint, let's recall what derivatives mean:

• First derivative $S'(x)$: This tells us the slope or rate of change of the function
• Second derivative $S''(x)$: This tells us how the slope is changing, or the "curvature" of the function

The Natural Cubic Spline Constraint:

The natural cubic spline constraint can be written as:

$$S''(t_0) = S''(t_{k+1}) = 0$$

where:

• $S''(x)$: the second derivative of the spline function with respect to $x$ (measures curvature)
• $t_0$: the first (leftmost) boundary knot (minimum value in the data range)
• $t_{k+1}$: the last (rightmost) boundary knot (maximum value in the data range)
• $k$: the number of interior knots (not including boundary knots)

Let's break this down:

• $S''(t_0) = 0$: The second derivative of the spline function at the left boundary equals zero (no curvature at the left edge)
• $S''(t_{k+1}) = 0$: The second derivative of the spline function at the right boundary equals zero (no curvature at the right edge)

What This Constraint Means:

When the second derivative equals zero, it means there's no curvature at that point. In other words, the function becomes linear (straight line) at the boundaries. This prevents the spline from "curving away" dramatically from the data, making extrapolation more stable and reasonable.

Benefits of Natural Cubic Splines:

1. Stable extrapolation: The spline won't predict wildly unrealistic values beyond the data range
2. Fewer parameters: The constraint reduces the number of free parameters we need to estimate
3. Smoother behavior: The linear behavior at boundaries often matches our expectations about how real-world relationships should behave

This visualization clearly demonstrates the practical importance of natural cubic splines when extrapolation is necessary. In the shaded extrapolation zones, the standard spline exhibits dramatic curvature that may produce unrealistic predictions, while the natural spline maintains linear behavior due to the zero second derivative constraint at the boundaries. Both methods fit the interior training data similarly, but their boundary conditions lead to vastly different extrapolation behavior. This makes natural cubic splines the preferred choice for applications requiring predictions beyond the training data range, such as forecasting or physical modeling where linear extrapolation is more defensible than arbitrary polynomial behavior.

The spline regression model maintains several important mathematical properties that make it both theoretically sound and computationally efficient.

The Least Squares Solution:

Just like in regular linear regression, we can find the optimal coefficients using the normal equation:

$$\hat{\boldsymbol{\beta}} = (\mathbf{B}^T\mathbf{B})^{-1}\mathbf{B}^T\mathbf{y}$$

where:

• $\hat{\boldsymbol{\beta}}$: the estimated coefficient vector of dimension $p \times 1$ (optimal weights for basis functions)
• $\mathbf{B}$: the spline basis matrix of dimension $n \times p$ (evaluated basis functions)
• $\mathbf{B}^T$: the transpose of $\mathbf{B}$ with dimension $p \times n$ (rows and columns swapped)
• $\mathbf{B}^T\mathbf{B}$: a $p \times p$ matrix representing basis function overlaps (Gram matrix)
• $(\mathbf{B}^T\mathbf{B})^{-1}$: the inverse of $\mathbf{B}^T\mathbf{B}$ with dimension $p \times p$
• $\mathbf{B}^T\mathbf{y}$: a $p \times 1$ vector representing projection of responses onto basis functions
• $\mathbf{y}$: the response vector of dimension $n \times 1$ (observed values)

Key Properties of B-spline Basis Functions:

The B-spline basis functions have several mathematical properties that make them particularly well-suited for regression:

• Compact support: Each basis function $B_j(x)$ is non-zero only over a small interval, typically spanning $d+1$ knot intervals (where $d$ is the degree). This means that for any given $x$, only a few basis functions contribute to the final function value, making computations efficient.

• Partition of unity: The basis functions sum to 1 at every point: $\sum_{j=1}^{p} B_j(x) = 1$ for all $x$ in the domain. This ensures that our spline function behaves predictably and doesn't "blow up" to infinity.

• Non-negativity: Each basis function satisfies $B_j(x) \geq 0$ for all $x$ and $j$. This property helps ensure numerical stability and prevents oscillations in the fitted function.

• Continuity: The basis functions automatically ensure the desired level of continuity at the knots. For example, cubic B-splines ensure that the function, its first derivative, and second derivative are all continuous, creating smooth transitions between polynomial pieces.

Computational Advantages:

These properties make B-splines particularly well-suited for regression problems because they:

1. Reduce computational cost: The compact support property means we only need to evaluate a few basis functions for each data point
2. Ensure numerical stability: The non-negativity and partition of unity properties prevent numerical issues
3. Provide smooth fits: The continuity properties ensure that the resulting function is smooth and well-behaved
4. Enable efficient matrix operations: The structure of B-splines makes the matrix $\mathbf{B}^T\mathbf{B}$ well-conditioned and easy to invert