Constructing Natural Cubic Splines

Why impose natural boundary constraints?

A cubic regression spline can bend flexibly between its knots. But a cubic polynomial continued toward the edge of the data range can behave in ways that are driven more by its mathematical form than by the observations. A natural cubic spline addresses the boundary by requiring zero second derivative at each boundary knot. This tutorial develops that definition and shows exactly how the two constraints reduce the number of free basis functions.

Let \(a\) and \(b\) be the boundary knots, with \(a<b\), and let \(\xi_1,\ldots,\xi_K\) be distinct interior knots satisfying \(a<\xi_1<\cdots<\xi_K<b\). We work with cubic splines that are twice continuously differentiable at the interior knots, as in the earlier material on cubic regression splines and knot continuity.

For regression, the mean response is represented as \(E(Y_i\mid x_i)=f(x_i)\). The spline construction concerns the form of \(f\), rather than a particular error distribution. When fitting by ordinary least squares, the usual regression model writes \(Y_i=f(x_i)+\varepsilon_i\); for classical standard errors, one commonly assumes independent errors with mean zero and constant variance, and normal errors when exact small-sample \(t\)- or \(F\)-inference is required. None of those error assumptions changes the boundary constraints or the basis derivation.

Starting from the cubic truncated-power basis

With \(K\) interior knots, a general cubic spline can be written in the truncated-power basis. Define \((u)_+=\max(u,0)\), so that \((x-\xi_j)_+^3\) is zero to the left of knot \(\xi_j\) and equals \((x-\xi_j)^3\) to its right. The general representation is

There are \(K+4\) coefficients: four polynomial coefficients and one coefficient for each interior knot. Its second derivative is

At the lower boundary \(a\), every interior-knot term is zero because \(a<\xi_j\). At the upper boundary \(b\), each truncated linear term is \(b-\xi_j\). Therefore the natural constraints become two linear equations in the coefficients:

How the two constraints reduce the basis

The first equation gives \(\beta_2=-3a\beta_3\). Substituting this into the second equation and solving for \(\beta_3\) gives

Because \(b-a\ne 0\), these equations determine \(\beta_2\) and \(\beta_3\) from the knot coefficients. The coefficients \(\beta_0,\beta_1\) and \(\gamma_1,\ldots,\gamma_K\) remain free. Thus the natural spline has \(K+2\) free coefficients, two fewer than the unrestricted cubic spline.

Substituting the solved coefficient into the original function gives a reduced basis representation. For each interior knot, define

Then every natural cubic spline with these knots can be written as

The reduced set of basis functions is therefore \(1,x,g_1(x),\ldots,g_K(x)\). The two ordinary polynomial columns \(x^2\) and \(x^3\) no longer appear as independent columns: their coefficients have been determined by the boundary constraints. The resulting basis has \(K+2\) columns, matching the coefficient count.

Worked constructions

Using the basis in regression

Once the knots and boundary interval have been chosen, the reduced basis can be evaluated at each observed predictor value. The regression design matrix has one column for the constant function, one for \(x\), and one for each \(g_j(x)\). Ordinary least squares estimates the \(K+2\) coefficients in that representation, provided the design matrix has full column rank. This is a standard linear regression in its coefficients even though the fitted mean can be curved as a function of \(x\).

The boundary knots are part of the specification, not merely bookkeeping. The constraints apply at those chosen values \(a\) and \(b\); changing them changes the basis functions and the spline space. In practice they are commonly set at the ends of the predictor range being modeled. The interior knots control where the cubic pieces join, while the natural constraints reduce the degrees of freedom available for the curve.

Common mistake

A frequent error is to read “natural boundary condition” as \(f(a)=f(b)=0\), or to force the fitted curve through specified endpoint responses. Those are constraints on function values, not on curvature. The natural cubic spline conditions are specifically \(f''(a)=0\) and \(f''(b)=0\). In the truncated-power derivation, it is the second-derivative equations that determine \(\beta_2\) and \(\beta_3\); \(\beta_0\) and \(\beta_1\) remain free to be estimated from the data.

Another mistake is to subtract two coefficients from the model without showing why. The reduction is justified because the two boundary equations are independent when \(a<b\), and the derivation explicitly solves them. Counting alone is not a substitute for applying the constraints to the chosen spline representation.