Splines & Smoothing Splines
Piecewise polynomials, knots, B-splines, and the roughness-penalty trade-off.
On this page
Beginner: the previous lesson's problem was a single polynomial trying to cover the entire x-range at once — push its degree up to capture some local wiggle in the middle of the data, and the whole curve, including the parts far from that wiggle, starts swinging wildly. A spline fixes this with a simple change of strategy: stop asking one curve to do everything. Chop the x-axis into pieces at a handful of chosen points called knots, fit a separate, low-degree (usually cubic) polynomial to each piece, and require the pieces to join up smoothly — matching value, slope, and curvature exactly at every knot, so the result looks like one continuous curve instead of a jagged patchwork of unrelated segments. Flexibility becomes local: a wiggle in one region no longer forces distortion everywhere else.
Intermediate: the version above — a small, hand-picked number of knots — is called a regression spline. A smoothing spline takes the idea further: instead of carefully choosing a handful of knot locations, put a knot at every single data point (as locally flexible as the data will ever allow) and then add a roughness penalty that punishes excessive curvature. One tuning knob, λ, now controls everything: crank it way up and the penalty crushes curvature to nothing, forcing the fit into a straight line; push it toward zero and the penalty disappears, letting the curve bend enough to pass through (interpolate) every point exactly. Between those two extremes sits a continuum of increasingly smooth fits, reached by turning one dial instead of re-choosing knot positions by hand.
Advanced: this is not a new design principle — it's the same loss + penalty template this module has used since ridge regression (2.3.2), just pointed at a different kind of object. Ridge and lasso penalize the size of a finite coefficient vector; a smoothing spline penalizes the curvature of an entire function. The loss term — sum of squared residuals — doesn't change at all.
Worked example. Imagine a week of noisy hourly temperature-sensor readings you want to denoise before further analysis. A regression spline with, say, 4 knots placed at sunrise/midday/sunset/midnight would need those positions chosen by hand (and chosen badly if the diurnal pattern shifts during the week). A smoothing spline sidesteps that choice entirely — it places a knot at every single hourly reading and lets λ alone decide how much of the raw jitter survives into the fitted curve, tuned once (often by cross-validation) rather than re-guessed by hand every time the data changes shape.
The first term is the familiar sum-of-squared-residuals loss. The second is the roughness penalty: f''(x) is the function's curvature at x, so squaring and integrating it over the whole domain gives a single number measuring total curvature — zero for a perfectly straight line, large for a wildly bending curve. λ trades the two off exactly like the module's general template, L(y, X\theta) + \lambda R(\theta) (see the Regression module overview) — except here the object being penalized, f, is an entire function rather than a coefficient vector θ, and R is measured by an integral of curvature rather than a vector norm.
Part 1 — the degrees-of-freedom count that forces degree 3. Suppose the spline has knots at n ordered points. Between consecutive knots there are n − 1 segments, and each segment gets its own cubic polynomial — 4 free coefficients (of 1, x, x², x³) per segment — for:
At each of the n − 2 interior knots (the two outer knots have nothing to match up beyond them), continuity is enforced three separate ways: the value from the left piece must equal the value from the right piece, their first derivatives (slopes) must match, and their second derivatives (curvatures) must match too. That's 3 constraints per interior knot:
Subtracting tells you how much freedom survives this bookkeeping:
That number, n + 2, is suspiciously close to n — and it should be. There are exactly n data values the spline still has to match, one per knot (the pieces on either side of an interior knot already agree there thanks to the continuity constraint above, so "match the y-value at knot i" is one new constraint, not two — even at the two boundary knots, only one piece touches each, so it's one constraint there too). Spending n of the n + 2 remaining degrees of freedom on those value-matches leaves exactly 2 left over — and a natural cubic spline spends those final two by declaring the second derivative to be exactly zero at the two outer boundary knots. Every unknown the algebra started with now has exactly one job, with nothing left over and nothing missing.
It's worth checking that cubic really is special here, not an arbitrary choice. Repeat the same count for degree-2 (quadratic) pieces: 3 coefficients per segment gives 3(n-1) unknowns, but continuity of value, slope, and curvature at each interior knot is still 3 constraints per knot (matching two full derivatives doesn't get cheaper just because the pieces themselves are simpler) — so the leftover freedom is 3(n-1) - 3(n-2) = 3, a constant that never grows with n. That's nowhere near enough left to also match n different data values once n grows past a handful — quadratic pieces glued this tightly collapse toward a single global quadratic with almost no local flexibility, defeating the entire point of using pieces. Degree 3 is the smallest degree where second-derivative continuity is affordable and enough freedom survives to actually fit the data — exactly why "cubic spline," not quadratic or quartic, is the name that stuck.
Part 2 — restricting to piecewise-cubic functions isn't an approximation. Every worked example so far has assumed the fitted function is piecewise-cubic from the start, as if that were a convenient simplification. It is not a simplification at all. A classical result in nonparametric regression (Reinsch, 1967) says something much stronger: among every twice-differentiable function f — polynomials, trigonometric functions, anything smooth enough for the penalty integral to make sense — the exact minimizer of the objective above is always a natural cubic spline with a knot at each of the n data points, for every value of λ ≥ 0.
The intuition, without the full calculus-of-variations proof: take any smooth candidate f, and let g be the natural cubic spline that interpolates f's own values at the n knots — a well-defined object regardless of what f looks like elsewhere. Because g agrees with f at every x_i by construction, the loss term is identical for f and g — swapping one for the other costs nothing on the data-fit side. On the penalty side, a classical fact about interpolation (the natural cubic spline interpolant has the least total curvature among every function passing through those same points) guarantees:
with equality only if f was already that spline. So g matches f's loss exactly and never has a larger penalty — g is at least as good as f, for every candidate f you could name. The true unconstrained optimum can therefore never be anything other than some natural cubic spline with knots at the data.
Where this is used: this is why every practical smoothing-spline implementation gets to sidestep what looks like an intractable optimization over an infinite-dimensional space of functions. Because the true minimizer is guaranteed in advance to live inside one specific, finite (n-dimensional) space — natural cubic splines with knots at the n data points — fitting one reduces to picking a basis for that finite space, writing the objective as a finite quadratic form in the basis coefficients, and solving a linear system. The function-space optimization and the finite-dimensional linear-algebra problem provably give the same answer, so software never has to solve the harder-looking one.
Blue: a smoothing spline, blended continuously between its two provable extremes — the natural-cubic-spline interpolant (λ→0) and the OLS line (λ→∞) — as the penalty-strength slider moves. Amber, dashed: a regression spline through only 2 fixed knots, fit by plain least squares with no penalty at all — drag its knots to reshape it, and compare how differently the two methods reach a similarly flexible-looking fit.
The same noisy data, refit as a plain interpolating natural cubic spline through only the K knots the slider chooses (dark circles); each independent cubic segment between two knots is drawn in its own color. The dashed tick at every interior knot is one straight line spanning both neighboring segments — visual proof that value and slope, not just the overall look, genuinely match exactly where the pieces meet, exactly the constraint the degrees-of-freedom count in the Derivation above spends its budget on.
Both from-scratch tabs build the truncated-power basis by hand and fit it with a hand-rolled Gaussian-elimination solver — a regression spline, exactly the "few fixed knots, no penalty" method from the diagram's amber curve. The library tab shows the real-world split this lesson keeps coming back to: UnivariateSpline's s parameter is a genuine roughness-penalty smoothing spline (a knot at every point), while SplineTransformer + LinearRegression is a regression spline (a handful of knots, ordinary least squares) — two different methods, both called "spline," living side by side in the same library.
#include <cmath>
#include <cstdio>
#include <vector>
#include <random>
#include <algorithm>
std::vector<double> basisRow(double x, const std::vector<double>& knots) {
std::vector<double> row = {1.0, x, x * x, x * x * x};
for (double k : knots) {
double d = x - k;
row.push_back(d > 0.0 ? d * d * d : 0.0);
}
return row;
}
// Solve A x = b with plain Gaussian elimination and partial pivoting --
// no linear-algebra library involved.
std::vector<double> solveLinearSystem(std::vector<std::vector<double>> A, std::vector<double> b) {
int n = static_cast<int>(b.size());
for (int col = 0; col < n; ++col) {
int pivot = col;
for (int r = col + 1; r < n; ++r)
if (std::fabs(A[r][col]) > std::fabs(A[pivot][col])) pivot = r;
std::swap(A[col], A[pivot]);
std::swap(b[col], b[pivot]);
for (int r = col + 1; r < n; ++r) {
double factor = A[r][col] / A[col][col];
for (int c = col; c < n; ++c) A[r][c] -= factor * A[col][c];
b[r] -= factor * b[col];
}
}
std::vector<double> x(n, 0.0);
for (int row = n - 1; row >= 0; --row) {
double total = b[row];
for (int c = row + 1; c < n; ++c) total -= A[row][c] * x[c];
x[row] = total / A[row][row];
}
return x;
}
int main() {
const int n = 40;
std::vector<double> knots = {2.5, 5.0, 7.5};
const int p = 4 + static_cast<int>(knots.size());
std::mt19937 rng(0);
std::normal_distribution<double> noise(0.0, 0.4);
std::vector<double> xs(n), ys(n);
for (int i = 0; i < n; ++i) {
xs[i] = i * 10.0 / (n - 1);
ys[i] = 4.0 + 2.2 * std::sin(xs[i] / 1.8) + noise(rng);
}
std::vector<std::vector<double>> XtX(p, std::vector<double>(p, 0.0));
std::vector<double> Xty(p, 0.0);
for (int i = 0; i < n; ++i) {
std::vector<double> row = basisRow(xs[i], knots);
for (int a = 0; a < p; ++a) {
Xty[a] += row[a] * ys[i];
for (int c = 0; c < p; ++c) XtX[a][c] += row[a] * row[c];
}
}
std::vector<double> theta = solveLinearSystem(XtX, Xty);
std::printf("Fitted coefficients:\n");
for (double t : theta) std::printf(" %.4f\n", t);
std::vector<double> queries = {1.0, 5.0, 9.0};
for (double xq : queries) {
std::vector<double> row = basisRow(xq, knots);
double yhat = 0.0;
for (int a = 0; a < p; ++a) yhat += theta[a] * row[a];
std::printf("x=%.2f -> fitted y = %.3f\n", xq, yhat);
}
return 0;
}- Smoothing noisy economic and sensor time series. Monthly unemployment figures, daily sensor readings, or any series with real signal buried in measurement noise gets a smoothing-spline pass to recover the underlying trend shape without forcing a hand-picked polynomial degree or knot count onto data whose wiggliness isn't known in advance.
- Growth curves in biology and medicine. Pediatric height/weight-for-age percentile charts (the kind used at checkups) are literally built by fitting smoothing splines to large reference populations — flexible enough to track the real, non-linear shape of childhood growth, smooth enough not to chase individual measurement noise.
- Non-parametric baseline-hazard estimation in survival analysis. Modeling how risk changes over time without assuming a fixed parametric shape (exponential, Weibull, …) is a natural job for a smoothing spline; a later Survival Analysis lesson returns to this in detail.
- Yield-curve construction in finance. Interest rates are only observed at a handful of bond maturities; smoothing splines interpolate a full, smooth term structure across all maturities from those sparse observations, without the instability a single high-degree global polynomial would introduce.
- Computer-aided design and computer graphics. The word "spline" originally named a physical flexible strip of wood or metal that draftsmen pinned through fixed points to draw smooth curves by hand — modern CAD software and font/vector-graphics rendering still use the mathematical descendants of that idea (splines and B-splines) to represent smooth curves and surfaces.
- Any exploratory curve fit where the right polynomial degree is unknown. Whenever a relationship looks non-linear but its exact functional form isn't known ahead of time, a smoothing spline gives a flexible, well-behaved fit tuned by a single cross-validated
λ, avoiding the previous lesson's problem of manually guessing a global polynomial degree.
- Choosing knot count or placement blindly. For a regression spline, too few knots underfits any sharp local feature the fixed pieces simply can't bend around; too many knots concentrated in one region can still locally overfit, even though the method as a whole has no global penalty to stop it. Data-driven knot placement (quantiles of the x-distribution, or cross-validating knot count) is standard practice for exactly this reason.
- Confusing "smoothing spline" with "regression spline." They are genuinely different methods that both happen to use the word "spline": a smoothing spline puts a knot at every data point and controls flexibility with a roughness penalty; a regression spline uses a small, fixed, hand-chosen number of knots and no penalty at all, fit by plain least squares (exactly the method in the code tabs above). Mixing up which one a paper or library function means is an easy, common error.
- Trusting extrapolation beyond the knot range. A natural cubic spline is, by construction, forced to continue linearly beyond its outermost knots — its second derivative is pinned to zero at both boundaries, and the natural spline is defined to keep that zero curvature forever past the edge of the data. That's a straight-line artifact of the boundary condition, not a real prediction about what happens out there — trusting the fitted curve's apparent shape just past the edge of the data is a mistake the construction itself doesn't support.
Going deeper
The truncated-power basis used in the code tabs above — 1, x, x², x³, plus one (x−k)₊³ term per knot — is the simplest possible basis for the space of piecewise-cubic functions with matching value/slope/curvature at each knot, which is exactly why it's the clearest one to teach from. It is rarely what production spline-fitting code actually uses. B-splines represent that same function space — same set of achievable curves, same degrees of freedom — with a different, much better-behaved basis: each individual B-spline basis function is nonzero over only a few adjacent segments (compact support), rather than the truncated-power basis's columns like x³ and (x−k)₊³, which look increasingly similar to each other as more knots are added and become highly correlated. That correlation makes the normal-equations matrix XᵗX nearly singular — numerically ill-conditioned — while the B-spline version of the same matrix stays sparse and well-conditioned. This is a basis-choice trade-off, not a different modeling idea: same space of functions, same fitted curve in exact arithmetic, dramatically better numerical behavior in floating-point practice. It's exactly why sklearn.preprocessing.SplineTransformer in the code above builds a B-spline basis under the hood rather than a truncated-power one.
A natural cubic smoothing spline is fit to data spanning x from 0 to 10. What does the fitted curve look like at x = 15, and why? Separately, if the code tabs above swapped the truncated-power basis for a B-spline basis of the same degree and knots, would the fitted curve itself change?
At x = 15 the curve is exactly a straight line — the tangent line at the rightmost knot (x = 10), extended forever. This isn't approximate: a natural cubic spline's defining boundary condition pins its second derivative to zero at the two outermost knots, and by construction the spline is defined to continue with that same zero curvature (i.e. linearly) past the edge of the data. Swapping the basis, on the other hand, would change nothing about the fitted curve itself. Truncated-power and B-spline bases span the exact same finite-dimensional space of piecewise-cubic functions with matching value/slope/curvature at each knot, so a least-squares fit in either basis lands on the same curve in exact arithmetic — B-splines are chosen purely for much better numerical conditioning (each basis function has compact support, unlike the increasingly correlated columns of the truncated-power basis), not for representing a different or larger set of possible curves.
The previous lesson's global polynomial and this lesson's spline solve the same problem — fitting a flexible, non-linear curve — with opposite strategies for controlling wiggliness: one polynomial fights for global smoothness by capping its degree, a spline achieves it locally, piece by piece, with continuity constraints doing the work instead. A smoothing spline pushes that idea to its limit — a knot at every point, flexibility reined in by a roughness penalty instead of a knot count — and the classical result above guarantees that restricting to piecewise-cubic functions loses nothing: the true optimum of the penalized objective, searched over every conceivable smooth function, always turns out to live exactly there. The next lesson, Generalized Additive Models (2.3.8), builds directly on this: a GAM is a sum of smooth functions, one per feature, and "smooth function" there usually means exactly the spline machinery built in this lesson, fit feature-by-feature instead of once.