natural spline, unmkpp gives 'wrong' answer

I have the following data

x=[-3 -2 1 4 6]
y=[2 0 3 1 4]

And I want to calculate the coefficient of the corresponding natural(!) Spline which is given by

 $s(x_j) = a_0 + a_1*(x-x_j) + a_2*(x-x_j)^2 + a_3*(x-x_j)^3$

I obtain

 |    |     a_3 |     a_2 |     a_1 |    a_0 | intv    |
 | s0 |  0.5044 |       0 | -2.5044 | 2.0000 | [-3,-2] |
 | s1 | -0.2831 |  1.5132 | -0.9912 |      0 | [-2,1]  |
 | s2 |  0.2217 | -1.0351 |  0.4430 | 3.0000 | [1,4]   |
 | s3 | -0.1601 |  0.9605 |  0.2193 | 1.0000 | [4,6]   |

However when I am using matlab

clear all

 x=[-3 -2 1 4 6]
 y=[2 0 3 1 4]
 pspl=spline(x,[0 y 0])% a natural spline 
 [breaks,coefs,l,k,d] = unmkpp(pspl) 

I obtain

 |  2.0595 | -4.0595 |       0 | 2.0000 | 
 | -0.3796 |  2.1190 | -1.9405 |      0 |
 |  0.3003 | -1.2976 |  0.5238 | 3.0000 |
 | -0.5387 |  1.4048 |  0.8452 | 1.0000 |

I don't know how to interpret this, the only terms with coincide are that of the right column which must correspond to the terms $a_j$ but for the rest I don't know.

Any help is appreciated.

Thanks

Uwe Brauer

[Stephen23]

@Uwe Brauer: this time I formatted your code for you. In future you can do it yourself by selecting the code text and clicking the {} Code button.

[Uwe Brauer]

sorry and thanks that you reformatted the code. I will keep this in mind!

Answer 1

[John D'Errico]

Well, it gives the right answer, to a different problem than the one that you think you posed. :)

So first of all, if you read the help for spline:

"However, if Y contains two more values than X has entries, then the first and last value in Y are used as the endslopes for the cubic spline."

That is NOT the definition of a NATURAL spline, although I can see how one might mistake that easily enough. A natural cubic spline has SECOND derivatives zero at the ends.

Next, the coefficients that you see here:

pspl.coefs
ans =
       2.0595      -4.0595            0            2
     -0.37963        2.119      -1.9405            0
      0.30026      -1.2976      0.52381            3
     -0.53869       1.4048      0.84524            1

are listed with the constant term LAST, so in decreasing order. The first element in each row is the coefficient of the cubic term (as you show it, with the knot in that interval subtracted off).

So, if we differentiate the spline ONCE, then evaluate it at the end points, we see:

fnval(fnder(pspl,1),[-3 6])
ans =
            0  -6.6613e-16

The first derivatives are indeed zero at the ends, although not the second derivatives.

fnval(fnder(pspl,2),[-3 6])
ans =
       -8.119      -3.6548

If your goal is to generate a natural cubic spline, I'll need to think for a second. :) A smoothing spline will suffice, if we crank down on tol.

pspl = spaps(x,y,0)
pspl = 
  struct with fields:

        form: 'B-'
       knots: [-3 -3 -3 -3 -2 1 4 6 6 6 6]
       coefs: [2 1.1652 -2.174 6.1053 -0.99415 2.5731 4]
      number: 7
       order: 4
         dim: 1

It produces a result in B-spline form, but s we see, it is indeed a natural cubic spline.

fnval(fnder(pspl,2),[-3 6])
ans =
  -2.6645e-15  -4.4409e-16

I'm pretty sure there is a direct way to generate a natural cubic spline in pp form, but my mind is drawing a blank at the moment. Drat, I almost forgot about csape.

pp = csape(x,y,'var')
pp = 
  struct with fields:

        form: 'pp'
      breaks: [-3 -2 1 4 6]
       coefs: [4×4 double]
      pieces: 4
       order: 4
         dim: 1

fnval(fnder(pp,2),[-3 6])
ans =
     0     0

fnval(fnder(pp,1),[-3 6])
ans =
      -2.5044       2.1404

pp.coefs
ans =
      0.50439            0      -2.5044            2
     -0.28314       1.5132     -0.99123            0
      0.22173      -1.0351      0.44298            3
     -0.16009      0.96053       0.2193            1

So a cubic spline in pp form, with natural end conditions. Of course, if you were to use my SLM toolbox, a spline with natural end conditions is one of the many options.

[Uwe Brauer]

ok it seems that

pspl=csape(x, y, [2 2] )
[breaks,coefs,l,k,d] = unmkpp(pspl)

gives the correct answer.

[Uwe Brauer]

thanks for showing me the

fnval(fnder(pp,2),[-3 6])

command, very useful indeed. And you are right I should have read carefully about endslopes these are indeed not the second derivatives. sigh

BTW you talk about your SLM toolbox, forgive me my ignorance, but could you provide please a pointer? thanks

[John D'Errico]

NP. I should have included that before. And as I said, it is easy enough to mis-read, just assuming it would produce a natural spline.

https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

[Uwe Brauer]

Hi

thanks for your link, I will try tomorrow since on my laptop I only have a outdated 2010 version running.

I thought about it again, I think matlab documentation could be improved. Natural splines are the first to be mentioned in many textbook, so help spline should say a word about it, it should also provide a link to csape, it provides a link to splines but this is not enough. I will open a new thread on this, maybe somebody from mathworks, will notice.

Uwe Brauer

[John D'Errico]

Natural splines are also not that valuable. Yes, they get written about a lot, because they are easy to build, and fairly easy to understand.

Splines are what I call a metaphorical model. The idea is that a cubic spline is a good model for the deflection of a thin flexible beam, that passes through a set of points. If one presumes that the deflection of a thin, flexible beam forms a good approximation for the shape of your curve, even though the physics behind the shape of that curve has NOTHING to do with the shape of a true beam, if that is true, then a spline serves as a metaphor. Such a metaphor allows you to predict the behavior of your process, to the extent that the metaphor is a good approximation to the physical process.

Some metaphors serve as better approximations than others, and the quality of that approximation allows you to trust the predictions of the metaphorical model.

So how good are natural splines as a model for a general curve? Sometimes they are great, other times not as good. The example I like to use when I taught about the use of splines as models is to fit the curves sin(x) and cos(x) over the same domains, that is [0,2*pi]. How well will a natural cubic spline do on these two curves? The problem is a flexible beam uses what are called natural boundary conditions, coming from the variational calculus solution for the shape of that beam. This implies that the second derivatives at the end points are zero. That works great to fit sin(x) over the domain [0,2*pi], since the second derivatives at the ends of the domain are indeed zero. The problem is when you try to fit cos(x) over the same domain, you see the curve fits more poorly.

In real life, most curves that you will fit will not truly have zero curvature at the ends. So a natural cubic spline is actually not that great of a choice for some random curve that will be fit. In general, you will probably get better performance form other boundary conditions, and the BCs that seems to make the most sense for some random spline model are often the not-a-knot end conditions.

The point is that while students are taught about natural cubic splines as a jumping off point, it is actually not that great of a choice for your general spline. To be honest, it is actually a good idea that a natural cubic fit is somewhat hard to obtain, since if it was easy, then novices who want to use a spline fit will gravitate to what they learned about early in their studies. As well, for the novice who knows absolutely nothing, they will be told by their compatriots (who also know little beyond what they learned from others) they will be told to use a natural spline fit.

While you may think this does not happen it really does. Such memes propagate all over numerical methods. For example, people are constantly wanting to use tools like the determinant to test for singularity of a matrix or they want to use inv to solve linear systems of equations. Both of these things are bad things in terms of numerical analysis, but because inv and det are do easy to use, people assume they must be appropriate.

A long argument, but as I said, I don't think it is really that good of a thing to make it easy to get a natural cubic spline interpolant. It is in there, if you know where to look. And if you understand enough about splines to also know why it may arguably be a bad thing, then you can find it without too much effort. (As I showed above, there were two distinct ways I came up with to produce a natural cubic interpolant.)

[Uwe Brauer]

Thanks for this long and detailed answer, I agree to certain agree. But I also give an elementary course in numerical analysis (and so do many colleagues in my department). Due to time constraints we can only discuss natural splines, since as you say are the easiest to understand and construct.

So in my experience Matlab is used a lot in such an academic environment, with the idea to familiarize the students which enable them later in there professional or academic carrier to do something mathematical more serious. However all I am saying is this, for this sort of environment, the matlab documentation could be a bit improved and make easy access to easy to handle tools.

[Michael Reshko]

Why does "order" property in your cubic splines shows 4? Shouldn't cubic be order 3?

[John D'Errico]

Why? Whats in a name?

pp = spline(rand(1,5),rand(1,5))
pp = 
  struct with fields:

      form: 'pp'
    breaks: [0.16679 0.33499 0.46094 0.51324 0.66433]
     coefs: [4×4 double]
    pieces: 4
     order: 4
       dim: 1

Blame it on the designers of splines in MATLAB. See that the CUBIC spline created is shown to have order 4.

Really, a name can indicate anything you want it to indicate. Many people use the words order and degree to be synonyms. My choice was to be consistent with that chosen in MATLAB. I did not write spline, but I need to use it and its cousins.

https://en.wikipedia.org/wiki/Degree_of_a_polynomial

In there we even see it stated that degree of a univariate polynomial in x would essentially be the highest power of x. As well, it states that order was commonly used in the past as a synonym for degree.

However, here we see a statement about the order of a spline as:

https://en.wikipedia.org/wiki/Order_of_a_polynomial

"the order of a spline, either the degree+1 of the polynomials defining the spline or the number of knot points used to determine it."

And here:

https://en.wikipedia.org/wiki/Spline_(mathematics)

"If the n polynomials Pi each have degree at most m, then the spline is said to be of degree m (or is of order m+1)."

A name is just a name. But order and degree have been disambiguated, no longer twins.