I'm looking to better understand how a Octave built in function splinefit works. That itself is a wrapper around something on the MATLAB file exchange.

As I understand it at a high-level, this function will fit polynomials to segments of a signal (the boundaries of which are knots or breakpoints). I can do that easily using ordinary least squares. But this also makes sure that neighbouring fits are continuous and smooth (hence, the similarity to splines and the name "splinefit"). This function is distinct from the idea of spline interpolation in that not every data point is used as a knot, only a few, and in each segment the polynomial should be optimal in a least squares sense (subject to the constraints at the segments boundaries) . It's utility seems to be in avoiding overfitting and that it will give you coefficients allowing you to build a model that can deal with nonlinearities by being piecewise.

I don't have the best stats background so I've really struggled with knowing what to look for in my literature search on understanding this. I am specifically interested in it as applied to cubic polynomials.

I don't even know what the general method is most properly called. I think it might be an example of "piecewise linear regression".

I've found in texts, lecture notes and online examples of something like this which seems similar but not right.

In that approach, you are building a single model based on basis functions like

$f(X) = \sum_{m=1}^{M} \beta_m h_m(X)$

where for a cubic regression spline with two knots at $\xi_1$ and $\xi_2$ is something like







I managed to build up a matrix $H$ and solve it to estimate coefficients. It looks okay but this isn't the same. The splinefit package actually fully determines a cubic function for each segment - i.e. I get four coefficients for each segment, as opposed to getting 3+number of knots coefficients to describe the whole signal.

I'd really appreciate any insight people can give me, or some direction to specific articles or literature. I'm having trouble filtering out the "noise" from similar yet different ideas.

Thanks in advance

  • $\begingroup$ The Numerical Recipes discussion of cubic spline interpolation, although not exactly the same as using splines for regression, is particularly lucid (yet brief) and presents the underlying concepts clearly. Copies of this section are available by Googling "numerical recipes spline" $\endgroup$
    – whuber
    Aug 1, 2020 at 12:34

2 Answers 2



The Matlab code you're referencing is using an over-parametrization with constraints instead of an unconstrained parametrization. What @Frank Harrel said about "[methods that use 4 parameters per segment] would allow for discontinuities in the derivatives of the function and a zero-oder discontinuity" is only true for unconstrained parametrizations, so it does not apply to this code. You can still write all twice differentiable piecewise polynomials with four parameters on each interval if you require the parameters to satisfy certain additional constraints. Let's look at an example:

Example of Two Different Representations of a Piecewise Polynomial

Denote the indicator function of a set $T$ by $$ I_T(x) = \begin{cases} 1, & \text{ if } x \in T \\ 0, & \text{ if } x \notin T \\ \end{cases}, $$ and suppose we have two knots $\xi_1$ and $\xi_2$. Let $A = (- \infty, \xi_1]$, $B = (\xi_1, \xi_2]$, and $C = (\xi_2, \infty)$. Then we could define some functions $$ \begin{align} g_1 (x) = I_A \ , \ g_2 (x) = x I_A \ & , \ g_3 (x) = x^2 I_A \ , \ g_4 (x) = x^3 I_A \\ g_5 (x) = I_B \ , \ g_6 (x) = x I_B \ & , \ g_7 (x) = x^2 I_B \ , \ g_8 (x) = x^3 I_B \\ g_9 (x) = I_C \ , \ g_{10} (x) = x I_C \ & , \ g_{11} (x) = x^2 I_C \ , \ g_{12} (x) = x^3 I_C \\ \end{align} $$ Let's also denote the set of all piecewise cubic polynomials with break points at the points $\xi_1$ and $\xi_2$ by $\mathcal{P}$. This set $\mathcal{P}$ is in fact a 12 dimensional vector space with $\{ g_i \}_{i=1}^{12}$ as a basis. If you don't know what all of this means, it essentially means that we need $12$ parameters to describe any piecewise cubic polynomials with break points at the points $\xi_1$ and $\xi_2$, and it means that any such piecewise polynomial $p \in \mathcal{P}$ can be written as $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ for 12 parameters $\{ \beta_i \}_{i = 1}^{12}$. If we pick some arbitrary values of the $\beta_i$'s, we can plot one of these functions:

enter image description here

Notice that the function is not continuous (and hence not differentiable, or second differentiable), which makes sense because up to this point we've been talking about the space of all piecewise cubic polynomials with breakpoints $\xi_1$ and $\xi_2$ with no reference to continuity or differentiability. If we want to talk about only the space of continuous piecewise polynomials with breakpoints $\xi_1$ and $\xi_2$, then we can define $\mathcal{P}^0 := \mathcal{P} \cap \mathcal{C}^0 (\mathbb{R})$ where $\mathcal{C}^0 (\mathbb{R})$ means the space of all continuous functions on $\mathbb{R}$.

For a function $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ to be continuous, it would need to satisfy two constraints, namely that $$ p(\xi_1) = \lim_{x \to \xi_1^+} p(x) \quad \text{ and } \quad p(\xi_2) = \lim_{x \to \xi_2^+} p(x) $$ or in terms of the parameters $\beta_i$: $$ \begin{align} \beta_1 + \beta_2 \xi_1 + \beta_3 \xi_1^2 + \beta_4 \xi_1^3 & = \beta_5 + \beta_6 \xi_1 + \beta_7 \xi_1^2 + \beta_8 \xi_1^3 \quad \text{ and }\\ \beta_5 + \beta_6 \xi_2 + \beta_7 \xi_2^2 + \beta_8 \xi_2^3 & = \beta_9 + \beta_{10} \xi_2 + \beta_{11} \xi_2^2 + \beta_{12} \xi_2^3 \\ \end{align} $$

But $\mathcal{P}^0$ is not a 12-dimensional space like $\mathcal{P}$ is! Essentially, having to satisfy two constraints subtracts two from the dimension to make it a 10-dimensional space. Requiring differentiability would require two more constraints, making $\mathcal{P}^1 := \mathcal{P} \cap \mathcal{C}^1 (\mathbb{R})$ an 8-dimensional space, and requiring twice differentiability would require yet two more constraints, making $\mathcal{P}^2 := \mathcal{P} \cap \mathcal{C}^2 (\mathbb{R})$ a 6-dimensional space. That number should be familiar as the same number of basis functions $\{ h_i \}_{i=1}^6$ you gave in your question, because $\{ h_i \}_{i=1}^6$ is a basis precisely for the 6-dimensional space $\mathcal{P}^2$. Here are the precise constraints in terms of the parameters $\beta_i$ needed to represent $\mathcal{P}^2$:

\begin{align*} \beta_1 + \beta_2 \xi_1 + \beta_3 \xi_1^2 + \beta_4 \xi_1^3 & = \beta_5 + \beta_6 \xi_1 + \beta_7 \xi_1^2 + \beta_8 \xi_1^3 \\ \beta_5 + \beta_6 \xi_2 + \beta_7 \xi_2^2 + \beta_8 \xi_2^3 & = \beta_9 + \beta_{10} \xi_2 + \beta_{11} \xi_2^2 + \beta_{12} \xi_2^3 \\ \beta_2 + 2 \beta_3 \xi_1 + 3 \beta_4 \xi_1^2 & = \beta_6 + 2 \beta_7 \xi_1 + 3 \beta_8 \xi_1^2 \\ \beta_6 + 2 \beta_7 \xi_2 + 3 \beta_8 \xi_2^2 & = \beta_{10} + 2 \beta_{11} \xi_2 + 3 \beta_{12} \xi_2^2 \\ 2 \beta_3 + 6 \beta_4 \xi_1 & = 2 \beta_7 + 6 \beta_8 \xi_1 \\ 2 \beta_7 + 6 \beta_8 \xi_2 & = 2 \beta_{11} + 6 \beta_{12} \xi_2 . \end{align*}

The first two are the continuity constraints from before, the next two are the two differentiability constraints for the points $\xi_1$ and $\xi_2$, and the last two are the second differentiability constraints.

But now we have two different ways of representing functions $p(x) \in \mathcal{P}^2$: we can either write them in terms of the functions $g_i (x)$ as $p(x) = \sum_{i=1}^{12} \beta_i g_i (x)$ where the $\beta_i$ are forced to satisfy the six constraints above, or we can write $p(x) = \sum_{i=1}^{6} \alpha_i h_i (x)$ for some different parameters $\alpha_i$ are aren't constrained at all.

This is the crux of the matter. The Matlab code is doing something similar using the 12 functions $g_i$ and then requiring the parameters to satisfy certain constraints. The set of functions $\{ g_i \}_{i=1}^{12}$ is no longer a basis for this space, because they won't be linearly independent, but you can call this set of functions a spanning set, redundant basis, or over-complete basis (the latter two are not actually bases at all, since they are not linearly independent). This language which is used frequently when talking about frames in signal processing and in dictionary learning.

B-Splines and Numerics

In the actual code, you won't see any functions like $g_i$ or $h_i$. The basis $\{ h_i \}_{i=1}^6$ is called a truncated power basis, and it has some problems. These start getting somewhat technical, so see Carl de Boor, A Practical Guide to Splines Chapter VIII "Example: the truncated power basis can be bad" for details if you are very interested in that. The solution is to use B-splines instead, which are a different basis for the same space with better numerical properties. The code you're referencing uses B-spline in its implementation, so if you want to fully understand it you'll have to learn about them as well.


Regression spline describes the concept of taking a fixed set of knots $\{ \xi_i \}_{i=1}^n$ and then fitting them to data according to some criteria, often involving least squares (the definition here is not given on Wikipedia, so see Hastie, Tibshirani, and Friedman The Elements of Statistical Learning Section 5.2 for this definition). Spline for regression and spline fitting also describe the same concept. These aren't popularly used terms, as you've discovered. You might notice that none of these terms have their own Wikipedia article (all three terms redirect to smoothing spline, which has a small note at the end on about regression splines). Smoothing splines are a particular type of regression splines designed to solve the knot-choosing problem, described below.

You said in your question that "This function is distinct from the idea of spline interpolation in that not every data point is used as a knot, only a few." In practice we can choose as many knots as we want and we can set them to be whatever we want. However, too much freedom can be a bad thing when fitting statistical models because of the over-fitting problem. To avoid this, one solution is to limit the number of knots. However, with a small number of knots, the particular choice of which knots to use becomes increasingly important. We now have an additional problem of making this choice. Sometimes theory can help us decide, in the few cases where there is theoretical reason to set a knot in a particular point. However, in many applications there is no good way to choose the knots, and your results will differ depending on the choice of the knots.

Smoothing splines are a related technique that solves this problem, which you can read about on the Wikipedia page.


Chapter 5 of Hastie, Tibshirani, and Friedman's The Elements of Statistical Learning is an excellent and short source on this. Carl de Boor's A Practical Guide to Splines (if you can find a copy) is the source on splines. Most of the book deals more with splines for interpolation, but chapter XIV is dedicated to smoothing splines.

  • $\begingroup$ Thanks for taking the time to give this excelent answer @erikperkerson, it has really helped me understand this a lot more. $\endgroup$ Aug 19, 2020 at 16:43

There are some great references and interactive demonstrations of splines. I refer to some of them in my course notes. I don't know of a method that uses 4 parameters per segment, because this would allow for discontinuities in the derivatives of the function and a zero-order discontinuity where you would see a vertical jump in the function value. To constrain the function to "meet" from the left and the right of a knot, and to constrain the first and second derivatives to also meet (i.e. to make it so that human eye cannot see the join point) you add only the one term like $(X - a)_{+}^{3}$ for each knot crossing in a cubic spline function. My notes and book go further to push the use of a linear tail-restricted cubic spline to reduce the total number of parameters by 4 and to prevent unreasonable wiggling on the left and right of the function.

  • $\begingroup$ Thanks @frank . Your lecture notes helped me "roll my own" which seems to do as good a job so this has helped me a lot (I am trying to port an Octave program that used splinefit in a step to C++, so I have substituted that in). I've not accepted it as the answer, just in case someone else can tell me more specifically what is happening in the octave code itself (like you say, it seems strange its getting four parameters per segment with continuity - I'd like to know how it does that!) $\endgroup$ Aug 4, 2020 at 6:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.