# TPRS Thin Plate Regression Splines (TPRS). Usage ``` python 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](TPRS.md#whittaker.TPRS.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](TPRS.md#whittaker.TPRS.d) is the covariate dimension. Common choices: `m=2` for `d <= 3` (the default), `m=3` for [d](TPRS.md#whittaker.TPRS.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: 1. **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`). 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 ``` python 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](#d) | Covariate dimension inferred from the training data. | | [eigenvalues](#eigenvalues) | Eigenvalues of the projected thin plate spline kernel matrix. | | [is_fitted](#is_fitted) | `True` after `fit()` has been called. | | [n_basis](#n_basis) | Total number of basis functions `k`. | | [null_space_dim](#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](TPRS.md#whittaker.TPRS.d), [eigenvalues](TPRS.md#whittaker.TPRS.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()](#basis_matrix) | Evaluate the TPRS basis at `x`. | | [fit()](#fit) | Fit the TPRS to training data `x`. | | [identifiability_constraints()](#identifiability_constraints) | Return the sum-to-zero constraint row for the intercept. | | [null_space_dimension()](#null_space_dimension) | Return the dimension of the unpenalized polynomial null space. | | [penalty_matrix()](#penalty_matrix) | Return the k x k penalty matrix S. | ------------------------------------------------------------------------ ### basis_matrix() Evaluate the TPRS basis at `x`. Usage ``` python 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 ``` python 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 ``` python 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 ``` python 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](TPRS.md#whittaker.TPRS.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 ``` python 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)`.