Spline Regression Methods

In polynomial regression, we expand a single predictor $x$ into a polynomial basis $(1, x, x^2, \ldots, x^d)$ and fit a linear model in these transformed features. This approach has a fundamental limitation: polynomial bases are global. Each basis function spans the entire range of $x$, meaning that the fitted value at any point depends on all the data—including data far away from that point.

This global nature creates two severe problems:

1. Runge's phenomenon: High-degree polynomials oscillate wildly near the boundaries, producing poor fits despite low training error.
2. Numerical instability: The Vandermonde matrix $(x_i^j)$ becomes extremely ill-conditioned as degree increases, causing coefficient estimates to explode.

Splines solve both problems by replacing global polynomial bases with piecewise polynomial bases that have local support—each basis function is non-zero only over a finite interval. The resulting fits adapt locally to the data while maintaining smoothness across the entire domain.

Before diving into B-splines specifically, we need to understand what a spline is in general.

Definition (Spline Function):

A spline of degree $d$ on the interval $[a, b]$ with knots $a = \xi_0 < \xi_1 < \xi_2 < \cdots < \xi_K < \xi_{K+1} = b$ is a function $s: [a, b] \to \mathbb{R}$ such that:

1. On each subinterval $[\xi_j, \xi_{j+1})$, the function $s$ is a polynomial of degree at most $d$.
2. The function $s$ has continuous derivatives up to order $d-1$ at each interior knot $\xi_1, \ldots, \xi_K$.

The knots $\xi_1, \ldots, \xi_K$ are the points where the polynomial pieces join. The smoothness requirement ensures that, despite being piecewise, the spline appears smooth to the eye—no visible corners or jumps up to the $(d-1)$-th derivative.

Degrees of Freedom of a Spline Space:

How many parameters are needed to specify a spline of degree $d$ with $K$ interior knots?

Consider: We have $K+1$ subintervals, each carrying a polynomial of degree $d$, which has $d+1$ coefficients. This gives $(K+1)(d+1)$ free parameters. However, at each of the $K$ interior knots, we impose $d$ continuity conditions (continuity of the function itself and its first $d-1$ derivatives). Thus:

$$\dim(\text{Spline Space}) = (K+1)(d+1) - Kd = K + d + 1$$

This is a fundamental formula. For cubic splines ($d = 3$) with $K$ interior knots, the dimension is $K + 4$. This matches: we need to specify $K$ knot positions plus 4 boundary conditions (or equivalently, choose from a $(K+4)$-dimensional basis).

The Basis Function Approach:

Since the space of splines of degree $d$ with $K$ knots has dimension $K + d + 1$, we can represent any such spline as a linear combination of $K + d + 1$ basis functions:

$$s(x) = \sum_{j=1}^{K+d+1} \beta_j B_j(x)$$

The question is: what basis should we use? The most natural-seeming basis—truncated power functions—turns out to be poorly suited for computation. Enter B-splines.

The most straightforward basis for the space of degree-$d$ splines with knots $\xi_1, \ldots, \xi_K$ is the truncated power basis:

$${1, x, x^2, \ldots, x^d, (x - \xi_1)+^d, (x - \xi_2)+^d, \ldots, (x - \xi_K)_+^d}$$

where the truncated power function is defined as:

$$(x - \xi)_+^d = \begin{cases} (x - \xi)^d & \text{if } x > \xi \ 0 & \text{if } x \leq \xi \end{cases}$$

This basis is intuitive: the first $d+1$ functions are the standard polynomial basis, and each additional function 'kicks in' at its corresponding knot to allow a new polynomial piece.

Why the Ill-Conditioning?

The problem is that truncated power functions are nearly collinear. Consider two adjacent knots $\xi_j$ and $\xi_{j+1}$. The functions $(x - \xi_j)+^d$ and $(x - \xi{j+1})_+^d$ are highly correlated—both are zero below their respective knots, and both grow polynomially above them, differing only by a slight phase shift.

Mathematically, if the knots are closely spaced, the truncated power functions become nearly linearly dependent, causing the design matrix to be nearly singular. The closer the knots, the worse the conditioning.

The Solution: B-Splines

B-splines provide an alternative basis for the exact same spline space that has excellent numerical properties. Instead of functions that 'switch on' at knots, B-splines are bell-shaped functions with compact support—each is non-zero only over a finite interval spanning a few adjacent knots.

B-splines (Basis splines) are a numerically stable basis for the spline function space. Unlike the truncated power basis, B-splines have compact support and are defined via an elegant recursive formula.

Extended Knot Sequence:

To construct B-splines of degree $d$ with $K$ interior knots $\xi_1 < \cdots < \xi_K$ on $[a, b]$, we first define an extended knot sequence by adding $d+1$ copies of the boundary values:

$$\underbrace{t_0 = t_1 = \cdots = t_d = a}{d+1 \text{ copies}}, \quad t{d+1} = \xi_1, \quad t_{d+2} = \xi_2, \quad \ldots, \quad t_{d+K} = \xi_K, \quad \underbrace{t_{d+K+1} = \cdots = t_{2d+K+1} = b}_{d+1 \text{ copies}}$$

This extended sequence has $2d + K + 2$ knots in total. The repeated boundary knots ensure that the B-splines interpolate exactly at the endpoints.

The Cox-de Boor Recurrence Formula:

B-splines are defined recursively. For notational convenience, let $B_{j,k}(x)$ denote the $j$-th B-spline of degree $k$ (so $k$ goes from 0 to $d$).

Base case (degree 0): $$B_{j,0}(x) = \begin{cases} 1 & \text{if } t_j \leq x < t_{j+1} \ 0 & \text{otherwise} \end{cases}$$

This is simply the indicator function for the interval $[t_j, t_{j+1})$—a rectangular pulse.

Recursive case (degree $k \geq 1$): $$B_{j,k}(x) = \frac{x - t_j}{t_{j+k} - t_j} B_{j,k-1}(x) + \frac{t_{j+k+1} - x}{t_{j+k+1} - t_{j+1}} B_{j+1,k-1}(x)$$

Here, we adopt the convention that $0/0 = 0$ whenever a denominator vanishes (which happens when consecutive knots coincide).

This recurrence is sometimes also written as: $$B_{j,k}(x) = \omega_{j,k}(x) \cdot B_{j,k-1}(x) + \big(1 - \omega_{j+1,k}(x)\big) \cdot B_{j+1,k-1}(x)$$

where $\omega_{j,k}(x) = \frac{x - t_j}{t_{j+k} - t_j}$ is a local linear interpolation weight.

Understanding the Recurrence Geometrically:

The Cox-de Boor formula has a beautiful geometric interpretation:

1. Start with degree-0 B-splines: rectangular pulses on each knot interval
2. At each level, each B-spline is a weighted average of two degree-$(k-1)$ B-splines
3. The weights are linear functions that transition smoothly across the knot span
4. Each recursion step adds a degree of smoothness

By the time we reach degree $d$, each B-spline $B_{j,d}(x)$ is:

• Non-zero only on the interval $[t_j, t_{j+d+1}]$ (local support)
• A piecewise polynomial of degree $d$
• $C^{d-1}$ continuous everywhere (smooth)
• Non-negative everywhere
• Bell-shaped (single maximum)

B-splines possess a remarkable collection of properties that make them ideal for statistical modeling. These properties are not just mathematically elegant—they translate directly into computational and interpretive advantages.

Proof of Partition of Unity:

The partition of unity property is crucial and follows from induction on degree. For degree 0, the B-splines are indicator functions of disjoint intervals, so exactly one equals 1 at any point.

For degree $k$, suppose the property holds for degree $k-1$. Then: $$\sum_j B_{j,k}(x) = \sum_j \left[ \omega_{j,k}(x) B_{j,k-1}(x) + (1-\omega_{j+1,k}(x)) B_{j+1,k-1}(x) \right]$$

Re-indexing the second sum and using $\omega_{j,k}(x) + (1-\omega_{j,k}(x)) = 1$: $$= \sum_j B_{j,k-1}(x) \cdot \left[ \omega_{j,k}(x) + (1-\omega_{j,k}(x)) \right] = \sum_j B_{j,k-1}(x) = 1$$

by the induction hypothesis. The repeated boundary knots ensure this works at the domain boundaries.

For regression, we evaluate each B-spline basis function at each data point to form a design matrix. Given $n$ observations $x_1, \ldots, x_n$ and $p = K + d + 1$ B-spline basis functions, the design matrix $\mathbf{B}$ is:

$$\mathbf{B}{i,j} = B{j,d}(x_i) \quad \text{for } i = 1, \ldots, n \text{ and } j = 0, \ldots, p-1$$

The spline regression model is then: $$y_i = \sum_{j=0}^{p-1} \beta_j B_{j,d}(x_i) + \epsilon_i = \mathbf{B}_i \boldsymbol{\beta} + \epsilon_i$$

In matrix form: $\mathbf{y} = \mathbf{B} \boldsymbol{\beta} + \boldsymbol{\epsilon}$

This is simply linear regression with the B-spline design matrix! We can apply ordinary least squares, ridge regression, or any other linear model technique.

Using scipy.interpolate for B-splines:

In practice, you should use optimized library implementations rather than the recursive formula. In Python, the scipy.interpolate.BSpline class provides efficient evaluation using a variant of the de Boor algorithm. For design matrix construction specifically, patsy and statsmodels provide the bs() function analogous to R's splines::bs() function:

import patsy

# Create B-spline design with 5 interior knots, degree 3
B = patsy.dmatrix("bs(x, df=8, degree=3, include_intercept=True) - 1", 
                   data={"x": x})

The df parameter specifies the degrees of freedom (= number of basis functions), which equals $K + d + 1$ for open-boundary B-splines.

Let's work through a complete example of B-spline regression, from data generation to model fitting to interpretation. We'll compare B-spline fits with different numbers of knots to illustrate the flexibility-smoothness tradeoff.

Interpreting the Results:

Note that $R^2$ always increases with more knots (more parameters = lower training error), but Adjusted $R^2$ penalizes complexity. Beyond a certain point, adding knots hurts generalization.

The condition number remains small (under 100) even with 20 knots—B-splines maintain numerical stability where the truncated power basis would have failed catastrophically.

The naive recursive evaluation of B-splines via Cox-de Boor is inefficient—it recomputes the same lower-degree B-splines multiple times. de Boor's algorithm provides a much faster approach for evaluating a B-spline curve $s(x) = \sum_{j} \beta_j B_{j,d}(x)$ at a given point $x$.

The key insight is that at any point $x$, at most $d+1$ B-splines are non-zero (due to local support). De Boor's algorithm exploits this by building up the value iteratively using only the non-zero B-splines.

Algorithm Sketch:

1. Find the knot interval containing $x$: find $r$ such that $t_r \leq x < t_{r+1}$
2. Initialize: $d_j^{[0]} = \beta_{r-d+j}$ for $j = 0, \ldots, d$ (the $d+1$ relevant coefficients)
3. Iterate for $k = 1, \ldots, d$: $$d_j^{[k]} = (1 - \alpha_{j,k}) d_j^{[k-1]} + \alpha_{j,k} d_{j+1}^{[k-1]}$$ where $\alpha_{j,k} = \frac{x - t_{r-d+k+j}}{t_{r+1+j} - t_{r-d+k+j}}$
4. Result: $s(x) = d_0^{[d]}$

This computes $s(x)$ in $O(d^2)$ operations, compared to $O(2^d)$ for naive recursion.

B-splines are one of several ways to represent spline functions. Understanding the relationships between representations clarifies when each is appropriate.

Key Relationships:

1. Truncated Power ↔ B-spline: Same spline space, different basis. Any truncated power spline can be written as a B-spline and vice versa—just with different coefficients. B-splines are preferred computationally.

2. B-spline → Natural Spline: Natural cubic splines are B-splines with additional constraints (zero second derivative at endpoints). We'll cover these in detail on Page 2.

3. B-spline Derivatives: The derivative of a degree-$d$ B-spline is a degree-$(d-1)$ B-spline: $$\frac{d}{dx} B_{j,d}(x) = d \left( \frac{B_{j,d-1}(x)}{t_{j+d} - t_j} - \frac{B_{j+1,d-1}(x)}{t_{j+d+1} - t_{j+1}} \right)$$

This is important for penalized splines (Page 5).

4. B-spline Integrals: Integration of B-splines leads to higher-degree B-splines. This is useful for computing penalty matrices.

We've covered substantial ground in establishing B-splines as the foundation for spline regression. Let's consolidate the key ideas:

What's Next:

With B-splines as our foundation, we can now explore more specialized spline constructions:

• Page 2: Natural Splines — Constrained cubic splines that reduce degrees of freedom and improve behavior at boundaries
• Page 3: Smoothing Splines — An optimization-based approach that avoids explicit knot selection entirely
• Page 4: Knot Selection — Strategies for choosing the number and placement of knots
• Page 5: Penalized Splines — Combining dense knot placement with penalty terms for the ultimate in flexibility