Splines and Smoothing
A straight line is often the wrong shape for \(E[Y \mid X = x]\). This page collects the ways of letting the regression function \(f\) bend in one predictor while the fit stays an ordinary least-squares problem (or a ridge-like one): polynomials, step functions, regression splines, natural splines, penalised splines and smoothing splines. All of them are linear smoothers, \(\hat{\mathbf y} = \mathbf S \mathbf y\) with \(\mathbf S\) depending only on the \(x\)’s, and that single fact carries the degrees-of-freedom bookkeeping and the cross-validation shortcut. 𐃏
The local alternative, fitting a fresh weighted regression at each query point, is Locally Weighted Regression. Splines are the global-but-piecewise answer; they feed straight into Additive Models, where one smooth is fitted per predictor.
polynomials and step functions
The non-parametric model is \(y_i = f(x_i) + \varepsilon_i\) with zero-mean, uncorrelated errors of variance \(\sigma^2\), and \(f\) left unspecified. Everything below works by choosing fixed functions of \(x\) and regressing on them.
Given known functions \(b_1, \dots, b_K\), fit \[ y_i = \beta_0 + \beta_1 b_1(x_i) + \dots + \beta_K b_K(x_i) + \varepsilon_i \] by least squares, with design-matrix columns \(b_j(x_1), \dots, b_j(x_n)\). Polynomial regression takes \(b_j(x) = x^j\); step-function regression takes \(b_j(x) = I(c_j \le x < c_{j+1})\).
The model is linear in \(\boldsymbol\beta\), so \(\hat{\boldsymbol\beta} = (\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top \mathbf y\) and all the inference machinery carries over: \[ \operatorname{se}\hat f(x_0) = \sqrt{\mathbf x_0^\top \widehat{\operatorname{Var}}(\hat{\boldsymbol\beta})\,\mathbf x_0} \] and \(\hat f(x_0) \pm 2\operatorname{se}\) is a pointwise band.
polynomial regression
A degree-\(p\) polynomial is flexible for large \(p\), but it imposes a global structure (data at age 30 move the curve at age 70), its coefficients mean little (R’s poly() uses orthogonal polynomials, so they are not even the \(\beta_j\) of \(x^j\)), and its variance explodes at the edges of the \(x\) range, where a high-degree tail is fitted by a handful of points.
lm(wage ~ poly(age, 4)). The fit and its standard error:
| age | \(\hat f\) | se |
|---|---|---|
| 18 | 51.9 | 5.30 |
| 20 | 64.6 | 3.59 |
| 30 | 103.7 | 1.41 |
| 40 | 117.4 | 1.05 |
| 50 | 119.6 | 1.13 |
| 60 | 116.6 | 1.83 |
| 70 | 107.0 | 3.69 |
| 80 | 81.6 | 13.63 |
The curve peaks near 48. The standard error at 80 is thirteen times the one at 40: only 72 of the 3000 men are over 64.5. For the binary response \(I(\text{wage} > 250)\) the same quartic inside a logit gives \(\hat\pi(40) = 0.034\) (band \(0.024\) to \(0.048\)) but a band of \([0, 0.73]\) at age 80. Bands for a logistic fit are built on the logit scale and transformed, which keeps them inside \((0,1)\).
step functions
Cut the range at \(c_1 < \dots < c_K\), define indicators \(C_0(x) = I(x < c_1)\), \(C_j(x) = I(c_j \le x < c_{j+1})\), \(C_K(x) = I(c_K \le x)\), and regress on \(C_1, \dots, C_K\). \(C_0\) is dropped because \(\sum_j C_j \equiv 1\) duplicates the intercept (the dummy-variable trap).
With \(\bar y_j\) the mean response in region \(j\), the least-squares estimates are \(\hat\beta_0 = \bar y_0\) and \(\hat\beta_j = \bar y_j - \bar y_0\). The fitted function is the regional sample mean, piece by piece.
lm(wage ~ cut(age, 4)) cuts 18–80 into four bands of width 15.5:
| region | \(n_j\) | \(\hat\beta\) (se) | \(\bar y_j\) |
|---|---|---|---|
| \((17.9, 33.5]\) (base) | 750 | 94.158 (1.476) | 94.158 |
| \((33.5, 49]\) | 1399 | 24.053 (1.829) | 118.212 |
| \((49, 64.5]\) | 779 | 23.665 (2.068) | 117.823 |
| \((64.5, 80.1]\) | 72 | 7.641 (4.987) | 101.799 |
Check: \(94.158 + 24.053 = 118.211 = \bar y_2\). The weakness shows in the first band: wages climb from about 52 to 110 between 18 and 33, but the step is flat at 94. Unless the predictor has natural breakpoints, a piecewise constant misses within-region trends.
proofs
Every observation in region 0 has fitted value \(\beta_0\), every one in region \(j \ge 1\) has \(\beta_0 + \beta_j\). Split the RSS by region:
\begin{align} \frac{\partial \text{RSS}}{\partial \beta_j} &= -2\sum_{i \in j}(y_i - \beta_0 - \beta_j) = 0 &&\Rightarrow\ \beta_0 + \beta_j = \bar y_j \quad (j \ge 1), \\ \frac{\partial \text{RSS}}{\partial \beta_0} &= -2\sum_{i \in 0}(y_i - \beta_0) - 2\sum_{j \ge 1}\sum_{i \in j}(y_i - \beta_0 - \beta_j) = 0 . \end{align}
The double sum vanishes by the first line, leaving \(\beta_0 = \bar y_0\), hence \(\hat\beta_j = \bar y_j - \bar y_0\).
regression splines
Splines combine the two ideas: cut the range at knots and fit a low-degree polynomial on each piece, then tie the pieces together.
piecewise polynomials and the constraint count
An unconstrained piecewise cubic with one knot at \(c\) is two separate cubics, \((K+1)(d+1) = 8\) parameters. On Wage with \(c = 50\) the left cubic reaches 116.95 at the knot and the right one starts at 112.93: a jump of \(-4.02\) that nothing in the model forbids.
Requiring continuity at the knot leaves 7 parameters; adding continuity of the first and second derivatives leaves 5. In general each independent linear constraint at a knot removes one dimension.
Re-expanding the right piece about the knot shows what a constraint looks like in a basis: \(f(x) = P_1(x) + \sum_{j=0}^{3}\delta_j (x-c)_+^j\) with \((x-c)_+ = \max(x-c, 0)\), and continuity of \(f, f’, f’’\) is exactly \(\delta_0 = \delta_1 = \delta_2 = 0\). What survives is \(1, x, x^2, x^3, (x-c)_+^3\).
lm(wage ~ threshold * poly(age, 3)) with threshold = I(age > 50)= looks like a “continuous” piecewise cubic but spans exactly the same 8-dimensional space as two separate cubics, so it reproduces the \(-4.02\) jump. Likewise threshold * age is two separate lines (jump \(-11.21\) at 50), not a linear spline. Continuity is a linear constraint on the coefficient vector; it has to be built into the basis. The genuine linear spline is lm(wage ~ age + pmax(age - 50, 0)):
\[ \hat f(x) = 60.13 + 1.3165\,x - 2.3221\,(x - 50)_+ , \]
slope \(1.3165\) before 50 and \(-1.0056\) after, meeting at the knot.
the truncated power basis
\(h(x, \xi) = (x - \xi)_+^3\) has \(h, h’, h’’\) continuous at \(\xi\) and \(h’’’\) jumping from 0 to 6. Adding \(\beta h(x, \xi)\) to a cubic breaks only the third derivative.
A cubic spline with knots \(\xi_1, \dots, \xi_K\) is a piecewise cubic with continuous first and second derivatives at every knot. The cubic splines form a vector space of dimension \(K + 4\) with basis \[ 1,\ x,\ x^2,\ x^3,\ (x - \xi_1)_+^3,\ \dots,\ (x - \xi_K)_+^3 , \] so fitting one is least squares on these columns. More generally a degree-\(d\) spline (\(d - 1\) continuous derivatives at each knot) has \(K + d + 1\) degrees of freedom and basis \(1, x, \dots, x^d, (x - \xi_k)_+^d\).
lm(wage ~ bs(age, knots = 50)) uses R’s B-spline basis (4 columns plus intercept); the truncated power basis gives the same curve to \(2.4 \times 10^{-10}\):
\[ \hat f(x) = -141.96 + 15.893\,x - 0.32167\,x^2 + 0.0021639\,x^3 - 0.0037213\,(x - 50)_+^3 . \]
Right of 50 the cubic coefficient becomes \(-0.0015574\) and the lower coefficients adjust so that value, slope and curvature match: \(\hat f(50) = 119.02\) from both sides, \(\hat f(18) = 52.5\), \(\hat f(80) = 78.3\). The picture is a single smooth arch with no visible join.
The two bases span the same space but behave very differently numerically. The columns \(x^3, (x - \xi_1)_+^3, (x - \xi_2)_+^3, \dots\) all grow like \(x^3\) for large \(x\) and are nearly collinear; B-splines are each non-zero over only a few adjacent knot intervals and stay well conditioned. Truncated-power coefficients are readable (each \(\beta_{3+k}\) is the change in the cubic coefficient at \(\xi_k\)); B-spline coefficients are not.
Colombia World Fertility Survey household counts by single year of age (0 to 98; 1581 at age 0, single digits after 89). A cubic spline with knots 25, 50, 75 has \(K + d + 1 = 7\) parameters. In the truncated power basis \[ \hat f(x) = 1392.80 + 80.516\,x - 6.9041\,x^2 + 0.11550\,x^3 - 0.14248\,(x-25)_+^3 + 0.035681\,(x-50)_+^3 - 0.016684\,(x-75)_+^3 , \] while the B-spline coefficients are \(1392.797\), \(670.968\), \(-863.818\), \(-968.094\), \(-1392.250\), \(-1373.655\), \(-1393.566\). Completely different numbers, identical fitted values (five-number summary \(-0.769\), \(82.18\), \(380.0\), \(920.2\), \(1657.7\), mean 563.1, for both). The intercepts agree because every other column in both bases vanishes at age 0. Knots at 20, 40, 60, 80 (8 parameters) or at the quartiles 24.5, 49, 73.5 give nearly the same curve: with this much smooth structure, knot placement barely matters.
proofs
Two cubics \(P_1, P_2\) carry 8 free coefficients. The constraints \(P_1( c) = P_2( c)\), \(P_1’( c) = P_2’( c)\), \(P_1’’( c) = P_2’’( c)\) are linear in the differences of coefficients with rows \((1, c, c^2, c^3)\), \((0, 1, 2c, 3c^2)\), \((0, 0, 2, 6c)\), which are linearly independent. Three independent equations on \(\mathbb R^8\) leave dimension 5, one at a time: \(8 \to 7 \to 6 \to 5\).
Left of \(\xi\) every derivative is 0. Right of \(\xi\): \(h = (x-\xi)^3\), \(h’ = 3(x-\xi)^2\), \(h’’ = 6(x-\xi)\), \(h’’’ = 6\). At \(\xi^+\) the first three are \(0, 0, 0\), matching the left; the third derivative is 6.
Each basis combination is a cubic spline: on every inter-knot interval each term is a cubic, and at \(\xi_k\) only \((x - \xi_k)_+^3\) is non-smooth, with two continuous derivatives by the lemma.
Dimension: \(K + 1\) cubic pieces carry \(4(K+1)\) coefficients and the knots impose \(3K\) independent constraints, leaving \(4(K+1) - 3K = K + 4\).
Independence: a vanishing combination is zero on \(x < \xi_1\), so its polynomial part is zero; then on \((\xi_1, \xi_2)\) it equals \(\beta_4 (x - \xi_1)^3\), so \(\beta_4 = 0\); continue knot by knot. \(K + 4\) independent functions in a \((K+4)\)-dimensional space form a basis.
Degree \(d\): \((K+1)(d+1)\) coefficients minus \(d\) constraints per knot,
\begin{align} (K+1)(d+1) - Kd = K + d + 1 , \end{align}
and \((x - \xi)_+^d\) has \(d - 1\) continuous derivatives with a jump of \(d!\) in the \(d\)-th.
natural splines
Beyond the outer knots a cubic spline is fitted by the few most extreme points and has the same boundary-variance problem as a polynomial.
A cubic spline constrained to be linear below the smallest knot and above the largest.
A natural cubic spline with \(K\) knots (boundary knots included) has \(K\) degrees of freedom; with \(K\) interior knots plus two boundary knots it has \(K + 2\) coefficients including the intercept.
R’s bs() and ns() both omit the intercept (lm supplies it) and their df counts their own columns. For bs(x, df, degree = d), \(\text{df} = K + d\), so \(K = \text{df} - d\) interior knots at the quantiles \(\tfrac{1}{K+1}, \dots, \tfrac{K}{K+1}\) (default \(d = 3\)). For ns(x, df), \(\text{df} = K + 1\), so \(K = \text{df} - 1\) interior knots, with boundary knots at the range of \(x\). On Wage:
| call | interior knots | coefficients |
|---|---|---|
bs(age, df = 4) | 42 (median) | 5 |
ns(age, df = 4) | 33.75, 42, 51 | 5 |
bs(age, df = 6) | 33.75, 42, 51 | 7 |
bs(age, df = 6, degree = 2) | 32, 39, 46, 53 | 7 |
ns(age, df = 6) | 30, 37, 42, 48, 54 | 7 |
(A 5th-degree spline with 10 knots has \(10 + 5 + 1 = 16\) degrees of freedom; bs() would call it df = 15.)
lm(wage ~ bs(age, df = 4)) against lm(wage ~ ns(age, df = 4)):
| age | cubic | se | natural | se |
|---|---|---|---|---|
| 18 | 52.5 | 5.41 | 58.6 | 4.24 |
| 20 | 64.8 | 3.67 | 66.9 | 3.34 |
| 30 | 103.7 | 1.49 | 103.1 | 1.51 |
| 40 | 117.2 | 1.10 | 117.8 | 1.26 |
| 50 | 119.9 | 1.09 | 119.6 | 1.18 |
| 60 | 116.3 | 1.81 | 115.7 | 1.65 |
| 70 | 105.8 | 3.66 | 107.0 | 3.92 |
| 80 | 88.1 | 12.13 | 95.9 | 8.40 |
Indistinguishable from 25 to 65; at the ends the natural spline is a straight line, so it cannot curl down, and its band at 80 is 30% narrower. With 6 df and knots at 25, 40, 60 (cubic) against ns(age, df = 6) the effect is stronger: se at 80 is 16.8 against 9.9, and at 18 it is 9.5 against 5.2.
proofs
A cubic spline with \(K\) knots has \(K + 4\) coefficients. Linearity left of the first knot zeroes the quadratic and cubic coefficients of the leftmost piece (two constraints), and likewise on the right: \(K + 4 - 4 = K\). With \(K\) interior plus two boundary knots, \((K + 2) + 4 - 4 = K + 2\).
penalised regression splines
More knots buy flexibility at the price of variance. On Wage, 25 equally spaced knots (29 coefficients) cut the RSS by only 0.5% relative to one knot (\(4.745\times 10^6\) against \(4.769 \times 10^6\)) while the curve wobbles from knot to knot and swings to a negative wage at 80. Instead of choosing few knots, use many and shrink.
Take the linear spline \(f(x) = \beta_0 + \beta_1 x + \sum_{j=1}^K \beta_{1j}(x - \xi_j)_+\). Constraining the knot coefficients by \(\max_j |\beta_{1j}| < C\), \(\sum_j |\beta_{1j}| < C\) or \(\sum_j \beta_{1j}^2 < C\) smooths the fit; the \(L_1\) version is a lasso that deletes knots, the \(L_2\) version a ridge that shrinks them all.
Let \(\mathbf D = \operatorname{diag}(0, 0, 1, \dots, 1)\), the \((K+2)\times(K+2)\) matrix that is zero on \((\beta_0, \beta_1)\) and the identity on the knot coefficients, so \(\boldsymbol\beta^\top \mathbf D \boldsymbol\beta = \sum_j \beta_{1j}^2\). Minimising \(\lVert \mathbf y - \mathbf X\boldsymbol\beta\rVert^2\) subject to \(\boldsymbol\beta^\top\mathbf D\boldsymbol\beta \le C\) is equivalent to minimising \[ \lVert \mathbf y - \mathbf X \boldsymbol\beta \rVert^2 + \lambda\, \boldsymbol\beta^\top \mathbf D \boldsymbol\beta \] for some \(\lambda \ge 0\), with solution \[ \hat{\boldsymbol\beta}_\lambda = (\mathbf X^\top \mathbf X + \lambda \mathbf D)^{-1}\mathbf X^\top \mathbf y . \]
This is ridge regression with \(\mathbf I\) replaced by \(\mathbf D\): the intercept and slope go unpenalised because a straight line is not rough. At \(\lambda = 0\) it is the unpenalised \(K\)-knot spline; as \(\lambda \to \infty\) the knot coefficients vanish and the fit becomes the least-squares line. The number of knots then matters little, and \(\lambda\) (or, for unpenalised splines, the number and placement of knots) is chosen by cross-validation; on Wage the CV error curve flattens at about 3 df (natural) or 4 df (cubic).
proofs
\begin{align} Q(\boldsymbol\beta) &= \mathbf y^\top\mathbf y - 2\boldsymbol\beta^\top\mathbf X^\top\mathbf y + \boldsymbol\beta^\top\mathbf X^\top\mathbf X\boldsymbol\beta + \lambda\boldsymbol\beta^\top\mathbf D\boldsymbol\beta &&\text{expand} \\ \nabla Q &= -2\mathbf X^\top\mathbf y + 2(\mathbf X^\top\mathbf X + \lambda\mathbf D)\boldsymbol\beta = 0 &&\text{set gradient to zero} \\ \hat{\boldsymbol\beta}_\lambda &= (\mathbf X^\top\mathbf X + \lambda\mathbf D)^{-1}\mathbf X^\top\mathbf y . &&\text{p.d. plus p.s.d. is invertible} \end{align}
The Hessian \(2(\mathbf X^\top\mathbf X + \lambda\mathbf D)\) is positive definite, so this is the minimum. Lagrange: the constrained problem has Lagrangian \(\lVert\mathbf y - \mathbf X\boldsymbol\beta\rVert^2 + \lambda(\boldsymbol\beta^\top\mathbf D\boldsymbol\beta - C)\); the \(-\lambda C\) term does not move the minimiser, and each \(\lambda\) corresponds to \(C = \hat{\boldsymbol\beta}_\lambda^\top\mathbf D\hat{\boldsymbol\beta}_\lambda\).
smoothing splines
Penalised splines still ask us to place knots. Smoothing splines remove that choice by penalising roughness directly over all functions.
For \(\lambda \ge 0\), the smoothing spline is the minimiser over twice-differentiable \(g\) of \[ \sum_{i=1}^{N}\bigl(y_i - g(x_i)\bigr)^2 + \lambda \int g’’(t)^2\,dt . \]
\(g’’\) measures how fast the slope changes, so the penalty is small for nearly linear \(g\) and large for wiggly \(g\). At \(\lambda = 0\) any smooth interpolant is optimal (RSS \(= 0\)); as \(\lambda \to \infty\) the minimiser is forced to \(g’’ \equiv 0\), and among lines the least-squares line wins.
The smoothing spline is a natural cubic spline with knots at the distinct \(x_i\): piecewise cubic, continuous first and second derivatives at each knot, linear outside the extreme knots. It is not the least-squares fit in that basis but a shrunken version of it, with \(\lambda\) setting the shrinkage.
The infinite-dimensional problem therefore collapses to a finite one. Let \(b_1, \dots, b_N\) span the natural cubic splines with knots at the data and set \(\{\mathbf B\}_{ij} = b_j(x_i)\).
With \(g = \sum_j \beta_j b_j\) the criterion becomes \[ (\mathbf y - \mathbf B\boldsymbol\beta)^\top(\mathbf y - \mathbf B\boldsymbol\beta) + \lambda\,\boldsymbol\beta^\top\boldsymbol\Omega_B\boldsymbol\beta, \qquad \{\boldsymbol\Omega_B\}_{jk} = \int b_j’’(t)\,b_k’’(t)\,dt , \] solved by \(\hat{\boldsymbol\beta} = (\mathbf B^\top\mathbf B + \lambda\boldsymbol\Omega_B)^{-1}\mathbf B^\top\mathbf y\), with \(\hat g(x) = \sum_j b_j(x)\hat\beta_j\).
smoother matrix and effective degrees of freedom
With \(N\) knots the spline looks like it spends \(N\) parameters, but they are shrunk. The fitted values are \(\hat{\mathbf g} = \mathbf S_\lambda \mathbf y\) with \(\mathbf S_\lambda = \mathbf B(\mathbf B^\top\mathbf B + \lambda\boldsymbol\Omega_B)^{-1}\mathbf B^\top\), depending only on the \(x\)’s and \(\lambda\).
\(\mathrm{df}_\lambda = \operatorname{tr}(\mathbf S_\lambda)\). For a regression spline with \(M\) basis columns \(\mathbf B_\xi\), the hat matrix \(\mathbf H_\xi = \mathbf B_\xi(\mathbf B_\xi^\top\mathbf B_\xi)^{-1}\mathbf B_\xi^\top\) has \(\operatorname{tr}(\mathbf H_\xi) = M\), the parameter count; the trace generalises that count to shrinking smoothers.
Both \(\mathbf H_\xi\) and \(\mathbf S_\lambda\) are symmetric positive semidefinite. \(\mathbf H_\xi\) is idempotent with rank \(M\) (eigenvalues 0 or 1: it projects). \(\mathbf S_\lambda\) has rank \(N\) and \(\mathbf S_\lambda^2 \preceq \mathbf S_\lambda\): in the Reinsch form \(\mathbf S_\lambda = (\mathbf I + \lambda\mathbf K)^{-1}\), \(\mathbf K = \mathbf B^{-\top}\boldsymbol\Omega_B\mathbf B^{-1}\), its eigenvalues are \[ \rho_k(\lambda) = \frac{1}{1 + \lambda d_k} \in (0, 1] , \] where \(d_k \ge 0\) are the eigenvalues of \(\mathbf K\). It shrinks each eigen-direction of \(\mathbf y\).
\(\mathrm{df}_\lambda = \sum_{k=1}^N (1 + \lambda d_k)^{-1}\) decreases continuously from \(N\) at \(\lambda = 0\) (interpolation) to \(2\) as \(\lambda \to \infty\) (the least-squares line).
Compare the ridge formula on Regularised Regression, \(\mathrm{df}(\lambda) = \sum_j d_j^2/(d_j^2 + \lambda)\): the same shape of fractional shrinkage factors, with the eigenvalues of the roughness matrix playing the role of \(1/d_j^2\). The difference is the null space: ridge shrinks everything to zero, the smoothing spline leaves the two-dimensional space of lines untouched.
For any \(\lambda\), \[ \text{RSS}_{\text{cv}}(\lambda) = \sum_{i=1}^N \bigl(y_i - \hat g_\lambda^{(-i)}(x_i)\bigr)^2 = \sum_{i=1}^N\left(\frac{y_i - \hat g_\lambda(x_i)}{1 - \{\mathbf S_\lambda\}_{ii}}\right)^2 , \] so every leave-one-out fit comes from the single full fit. Generalised cross-validation replaces each \(\{\mathbf S_\lambda\}_{ii}\) by the average \(\mathrm{df}_\lambda / N\).
smooth.spline(age, wage, df = 16) solves \(\operatorname{tr}(\mathbf S_\lambda) = 16\), landing at \(\lambda = 6.5\times 10^{-4}\): a wiggly curve with small bumps at the age heaps. smooth.spline(age, wage, cv = TRUE) minimises the leave-one-out criterion, landing at \(\lambda = 0.028\) with \(\mathrm{df}_\lambda = 6.79\), nearly indistinguishable from ns(age, df = 6). R warns that CV “with non-unique x values seems doubtful”: the 3000 ages take only 61 distinct values, so the spline has \(N = 61\) knots and leaving out one \(y\) of a tie does not leave out its \(x\). The two fits differ little, so the 6.8-df spline is preferred.
proofs
Take any candidate \(g\) and let \(\tilde g\) be the natural cubic spline with knots at the \(x_i\) interpolating \(g(x_i)\) (it exists and is unique: that space has dimension \(N\)). Both have the same RSS. Put \(u = g - \tilde g\), which vanishes at every knot:
\begin{align} \int \tilde g’’\,u’’ &= \bigl[\tilde g’’ u’\bigr]_{x_1}^{x_N} - \int \tilde g’’’ u’ &&\text{by parts} \\ &= -\sum_{j=1}^{N-1}\tilde g’’’(x_j^+)\bigl[u(x_{j+1}) - u(x_j)\bigr] = 0 . &&\tilde g’’ = 0 \text{ at the ends; } \tilde g’’’ \text{ piecewise constant; } u = 0 \text{ at knots} \end{align}
So \(\int g’’^2 = \int \tilde g’’^2 + \int u’’^2 \ge \int \tilde g’’^2\): every candidate is beaten or matched by a natural cubic spline, and the minimiser is one. Within that family the criterion is the generalised ridge below.
\(g(x_i) = (\mathbf B\boldsymbol\beta)_i\) gives the loss. For the penalty,
\begin{align} \int g’’(t)^2\,dt = \int \Bigl(\sum_j \beta_j b_j’’\Bigr)\Bigl(\sum_k \beta_k b_k’’\Bigr)dt = \sum_{j,k}\beta_j\beta_k\int b_j’’ b_k’’\,dt = \boldsymbol\beta^\top\boldsymbol\Omega_B\boldsymbol\beta . \end{align}
\(\boldsymbol\Omega_B\) is symmetric p.s.d. (\(\boldsymbol\beta^\top\boldsymbol\Omega_B\boldsymbol\beta = \int g’’^2 \ge 0\)), so the penalised-spline derivation applies with \(\mathbf B, \boldsymbol\Omega_B\) for \(\mathbf X, \mathbf D\).
Let \(\mathbf A = \mathbf B^\top\mathbf B + \lambda\boldsymbol\Omega_B\) (p.d.). Symmetry is immediate; \(\mathbf v^\top\mathbf S_\lambda\mathbf v = \mathbf u^\top\mathbf A^{-1}\mathbf u \ge 0\) with \(\mathbf u = \mathbf B^\top\mathbf v\). \(\mathbf H_\xi^2 = \mathbf H_\xi\) by cancelling \((\mathbf B_\xi^\top\mathbf B_\xi)^{-1}\mathbf B_\xi^\top\mathbf B_\xi = \mathbf I\), and \(\operatorname{tr}\mathbf H_\xi = \operatorname{tr}\bigl((\mathbf B_\xi^\top\mathbf B_\xi)^{-1}\mathbf B_\xi^\top\mathbf B_\xi\bigr) = M\). With \(N\) knots \(\mathbf B\) is square and invertible, so
\begin{align} \mathbf S_\lambda = \bigl(\mathbf B^{-\top}(\mathbf B^\top\mathbf B + \lambda\boldsymbol\Omega_B)\mathbf B^{-1}\bigr)^{-1} = (\mathbf I + \lambda\mathbf K)^{-1}. \end{align}
Diagonalising \(\mathbf K = \mathbf U\operatorname{diag}(d_k)\mathbf U^\top\) gives eigenvalues \(\rho_k = 1/(1 + \lambda d_k) \in (0, 1]\); then \(\mathbf S_\lambda - \mathbf S_\lambda^2\) has eigenvalues \(\rho_k(1 - \rho_k) \ge 0\), and all \(\rho_k > 0\) gives rank \(N\).
The trace is the sum of the eigenvalues, each decreasing in \(\lambda\). At \(\lambda = 0\) every \(\rho_k = 1\). As \(\lambda \to \infty\), \(\rho_k \to 0\) when \(d_k > 0\) and stays 1 when \(d_k = 0\). The null space of \(\mathbf K\) is the set of \(\boldsymbol\beta\) with \(\int g’’^2 = 0\), i.e. \(g(x) = a + bx\), which has dimension 2. So \(\mathrm{df}_\lambda \to 2\).
Replace \(y_i\) by \(\tilde y_i = \hat g^{(-i)}(x_i)\). The leave-one-out spline fits this point with zero residual, so it is also the smoothing spline for the modified data: \(\hat g^{(-i)}(x_i) = (\mathbf S_\lambda\tilde{\mathbf y})_i\). Writing \(S_{ij}\) for the entries,
\begin{align} \hat g^{(-i)}(x_i) &= \sum_{j \ne i} S_{ij}y_j + S_{ii}\hat g^{(-i)}(x_i) = \hat g(x_i) - S_{ii}y_i + S_{ii}\hat g^{(-i)}(x_i) \\ \Rightarrow\ y_i - \hat g^{(-i)}(x_i) &= y_i - \frac{\hat g(x_i) - S_{ii}y_i}{1 - S_{ii}} = \frac{y_i - \hat g(x_i)}{1 - S_{ii}} . \end{align}
With \(\mathbf H\) in place of \(\mathbf S_\lambda\) this is the familiar LOOCV identity for least squares.
see also
- Locally Weighted Regression – the local alternative: a fresh weighted fit at every query point instead of a global piecewise basis
- Regularised Regression – penalised and smoothing splines are generalised ridge; the trace-of-smoother df is the same idea
- Additive Models – one spline per predictor, fitted by backfitting
- Model Selection – choosing knots, df or \(\lambda\) by cross-validation and information criteria
- bias-variance decomposition – what every knot, degree and \(\lambda\) here is trading
- Linear Regression – every regression spline is OLS on a cleverly chosen design matrix
Comments