TPRS

Thin Plate Regression Splines (TPRS).

Usage

TPRS(
    k=10,
    m=2,
)

A thin plate spline is the function that minimizes squared error subject to a penalty on total curvature, with no need to choose knot locations — it is the natural multivariate generalization of the cubic smoothing spline. The full thin plate spline has one basis function per unique data point, which is computationally impractical for anything beyond a few hundred observations, so TPRS instead constructs a low-rank approximation using the leading eigenvectors of the (nullspace-projected) thin plate spline kernel matrix (Wood 2003). Because it works for any number of covariate dimensions d and requires no knot placement, TPRS is a good default basis for smooth terms of one or more continuous, non-cyclic covariates, especially in more than two dimensions where tensor-product alternatives become unwieldy.

Parameters

k: int = 10: Total number of basis functions (including the M null-space columns). Must satisfy k > M. Larger k allows more wiggly fits at the cost of more computation; the penalty (not k) ultimately controls smoothness once lambda is chosen. The default is 10.
m: int = 2: Spline order. Controls the order of derivative penalized: m=2 penalizes (squared) second derivatives, the classic “thin plate” bending energy. Must satisfy 2m > d where d is the covariate dimension. Common choices: m=2 for d <= 3 (the default), m=3 for d in {4, 5}. The default is 2.

Notes

The full thin plate spline basis uses the radial kernel

\eta_m(r) = \begin{cases} r^{2m - d} & 2m - d \text{ odd} \\ r^{2m-d} \log(r) & 2m - d \text{ even} \end{cases},

evaluated at pairwise distances r = ||x_i - x_j|| between data points, plus a polynomial null space of all monomials of total degree at most m - 1, which has dimension M = C(m - 1 + d, d). TPRS builds the basis in two stages:

Polynomial null space (first M columns): the unpenalized low-degree polynomials, for which the roughness penalty is identically zero (e.g. any straight line has zero bending energy under m=2).
Truncated spline part (remaining k - M columns): the full kernel matrix E is projected onto the orthogonal complement of the polynomial null space and eigen-decomposed; the k - M eigenvectors with the largest eigenvalues give the best rank-(k - M) approximation to the full thin plate spline in the sense of minimizing the change in the penalty for a given basis dimension.

The resulting penalty matrix is block-diagonal,

\mathbf{S} = \operatorname{diag}(0, \ldots, 0, \lambda_1, \ldots, \lambda_{k-M}),

with the M null-space rows/columns exactly zero and the remaining diagonal entries equal to the retained eigenvalues of the projected kernel matrix. Because the basis is derived from an eigendecomposition of a matrix that mixes all covariate scales, columns of x with very different scales can cause numerical issues; centering and/or standardizing each column of x before fitting is advisable.

Examples

import numpy as np
from whittaker.smooths import TPRS

rng = np.random.default_rng(0)
x = rng.uniform(0, 1, 100)

basis = TPRS(k=10).fit(x)
B = basis.basis_matrix(x)
S = basis.penalty_matrix()
B.shape, S.shape

((100, 10), (10, 10))

Attributes

Name	Description
d	Covariate dimension inferred from the training data.
eigenvalues	Eigenvalues of the projected thin plate spline kernel matrix.
is_fitted	`True` after `fit()` has been called.
n_basis	Total number of basis functions `k`.
null_space_dim	Dimension of the unpenalized null space of the TPRS penalty.

d

Covariate dimension inferred from the training data.

d: int

Set during fit() from the number of columns of the (reshaped) training covariates x; used to validate that new data passed to basis_matrix() has a matching number of columns.

eigenvalues

Eigenvalues of the projected thin plate spline kernel matrix.

eigenvalues: NDArray

These are the k - M largest eigenvalues retained from the eigendecomposition performed in fit(); they populate the diagonal of the penalized block of penalty_matrix() and determine how strongly each retained spline basis function is penalized.

is_fitted

True after fit() has been called.

is_fitted: bool

basis_matrix(), penalty_matrix(), and the convenience properties (d, eigenvalues, etc.) all require the basis to have been fitted first and raise RuntimeError otherwise.

n_basis

Total number of basis functions k.

n_basis: int

Equal to the k argument supplied at construction, i.e. the combined size of the polynomial null space (M columns) and the truncated spline part (k - M columns).

null_space_dim

Dimension of the unpenalized null space of the TPRS penalty.

null_space_dim: int

This is a convenience property equivalent to calling null_space_dimension(). For a thin plate regression spline of order m on d covariates, the null space has dimension M = \binom{m + d - 1}{d}, corresponding to the polynomial terms of degree less than m that are left unpenalized.

Methods

Name	Description
basis_matrix()	Evaluate the TPRS basis at `x`.
fit()	Fit the TPRS to training data `x`.
identifiability_constraints()	Return the sum-to-zero constraint row for the intercept.
null_space_dimension()	Return the dimension of the unpenalized polynomial null space.
penalty_matrix()	Return the k x k penalty matrix S.

basis_matrix()

Evaluate the TPRS basis at x.

Usage

Source

basis_matrix(x)

Parameters

x: NDArray: Covariate values. Shape (n,) or (n, d) where d must match the training dimension.

Returns

NDArray: Design matrix of shape (n, k). The first M = d + 1 columns are the polynomial null-space functions; the remaining k - M columns are the truncated spline functions.

fit()

Fit the TPRS to training data x.

Usage

Source

fit(x)

Parameters

x: NDArray: Training covariates. Shape (n,) for univariate or (n, d) for multivariate.

Returns

TPRS: Returns self for method chaining.

Raises

ValueError: If k is too large for the number of observations or too small for the covariate dimension.

identifiability_constraints()

Return the sum-to-zero constraint row for the intercept.

Usage

Source

identifiability_constraints()

When a TPRS term is combined with other terms sharing an intercept, it must be constrained so its contribution has zero mean over the training data to remain identifiable. The constraint is C @ beta = 0 where C is the column mean of the training basis matrix, i.e. colMeans(B_train) @ beta = 0.

Returns

NDArray: Constraint matrix of shape (1, k).

null_space_dimension()

Return the dimension of the unpenalized polynomial null space.

Usage

Source

null_space_dimension()

The TPRS penalty is exactly zero on polynomials of total degree at most m - 1, which span an M-dimensional null space with M = C(m - 1 + d, d) (the number of monomials of degree <= m - 1 in d variables). These occupy the first M columns of the basis matrix.

Returns

int: The null-space dimension M.

penalty_matrix()

Return the k x k penalty matrix S.

Usage

Source

penalty_matrix()

S is block-diagonal:

First M rows/cols: all zeros (null-space, unpenalised).
Remaining k - M rows/cols: diagonal entries equal to the leading eigenvalues of the projected kernel matrix.

Returns

NDArray: Shape (k, k).