using polyval with one of the polynomials of a spline obtained with the unmkpp command, the results do not overlap! Why?

I first interpolated some data using the spline command. Then I calculated the coefficients of the different cubic polynomials that make up the spline. Then I used the polyval command to check the first polynomial but the spline and polyval curves do not overlap.

What am I doing wrong?.

There is definitely something related to the various polynomials that I have not understood correctly. I thank anyone who helps me understand better. See the file script please:

%butta_sto_test

%

clear all

clc

% Data
% Time
Time=[7, 9, 11, 12]
% Temperture
Temperature=[49, 57, 71, 75]
% Data Time query
Timeq=[7:0.1:12];

% Results
% piecewise polynomial structure
pp=spline(Time, Temperature);
% Timeq data interpolation
Temperatureq=spline(Time, Temperature, Timeq);
% Calculating the coefficients of the 3 polynomials that make up the spline
[breaks, coeffs, m, n] = unmkpp(spline(Time,Temperature));
% First polinomial p1
p1=coeffs(1,:)
% First time interval corresponding to the first polynomial of the spline
Timeq1=[7:0.1:9];
% using polyval to calculate the values corresponding to Timeq1 using the polynomial p1
Temperatureq1=polyval(p1,Timeq1);

% Spline plot
figure(1)
plot(Timeq, Temperatureq, Time, Temperature, 'o', Time, Temperature, '--')
hold on
plot(Timeq1, Temperatureq1,'*')
grid on

Answer 1

[David Goodmanson]

Hello TH,.

spline polynomials are local to the secgment in question starting from x = 0, so if an interval is, for example 3 to 4.4, polyval uses 0 to 1.4. If you use

Timeq1=[7:0.1:9];                   % existing line
Timeq1poly = Timeq1-Timeq1(1);
% using polyval to calculate the values corresponding to Timeq1 using the polynomial p1
Temperatureq1=polyval(p1,Timeq1poly);

then it works.

[Steven Lord]

That is correct. From the description of the pp output argument on the spline documentation page:

Since the polynomial coefficients in coefs are local coefficients for each interval, you must subtract the lower endpoint of the corresponding knot interval to use the coefficients in a conventional polynomial equation. In other words, for the coefficients [a,b,c,d] on the interval [x1,x2], the corresponding polynomial is

[T Hafid]

I understand perfectly. Thank you very much.

[T Hafid]

Thank you

Answer 2

[dpb]

% Data
% Time
Time=[7, 9, 11, 12];
% Temperture
Temperature=[49, 57, 71, 75];
% Data Time query
Timeq=[7:0.1:12];

% Results
% piecewise polynomial structure
pp=spline(Time, Temperature);
% Timeq data interpolation
Temperatureq=spline(Time, Temperature, Timeq);

plot(Time,Temperature,'*')
hold on
plot(Timeq,Temperatureq)
Ts=ppval(pp,Timeq);
plot(Timeq,Ts,'x')
[bk,p]=unmkpp(pp)

bk = 1x4

7 9 11 12

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

p = 3x4

-0.3500 2.8500 -0.3000 49.0000 -0.3500 0.7500 6.9000 57.0000 -0.3500 -1.3500 5.7000 71.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

T1=bk(1):0.1:bk(2);
Tp=polyval(p(1,:),T1-bk(1));
plot(T1,Tp,'+')
xlim([7 9]), ylim([45 60])

Spline coefficients are local to the break region so to evaluate by conventional polynomial functions you must subract the leading breakpoint from each section. Since this is done for you automagically by ppval, it's nonsense to not use the toolset for the job.

[David Goodmanson]

It may be nonsense to not use spline & ppval for a production-type of job, but it is the exact opposite of nonsense, it's commendable to have curiosity and make the comparison that the OP did here. If everything is treated like a black box provided by Mathworks (and they do not have that attitude), how do people gain knowledge and perspective about the mathematical processes they are using?

[dpb]

I didn't say it wasn't a worthwhile exercise to figure out; only to not use the wrong toolset for the job... :)

Perhaps "nonsense" wasn't a good choice, granted.

[T Hafid]

Thanks for the explanations

Answer 3

[John D'Errico]

@David Goodmanson and @dpb and @Steven Lord have all said the polynomial segments are intended to be LOCAL with respect to the break points. I'm adding this as an answer only because I want to explain why it is done that way.

Supose I gave you two curves to interpolate using a spline, thus f(x1,y) and f(x2,y).

x1 = 0:1:10;
x2 = 1e8 + (0:1:10);

y = rand(1,11); % I dont give a hoot what values you choose for y.

spl1 = spline(x1,y);
spl2 = spline(x2,y);

The two curves should, in theory, be virtually identical. All I have done is horizontally translate the curve in x by some very large number.

Remember, these will be cubic splines. If the cubic segments for these splines were not local things, then in spl2, you would be forming a cubic polynomial where you are raising numbers on the order of 1e8 to the third power, and then adding and subtracting those values to other nubmers. Those polynomials would be unmanagable. You would expect to see massive numerical problems.

You don't want that. Even smaller shifts would cause serious problems. The solution is trivial. Just formulate each segment to be defined as a cubic polynomial, but where the LEFT end of each segment ALWAYS starts at zero, and has length given by diff(breaks). This resolves the numerical problems as much as possible.

[dpb]

@John D'Errico. I have wondered why the builtin tools don't also use the standardized centering and unit standard deviation as can do with polyfit/polyval, etc., ...

I suppose there are some implementatiions that do.

[John D'Errico]

@dpb That is an entirely valid question. I can offer a few reasons, some more valid than others.

1. The author did not think of it,
2. The author did not think it necessary.
3. Use of centering and scaling adds some slight additional mathematics to compute the polynomial coefficients. In my eyes, the choices would be to shift/scale things so the polynomial lives always on [-1,1], but I could also see arguments to map the interval always to [0,1]. I might like the latter, personally, but the difference is arguably minimal. A [0,1] interval makes the polynomial a bit easier to inspect, just by looking at the coefficients. And since that is something I have actually had reason to do on occasion, I might instinctively prefer it. But, do we need it? See below.

Alternative 1 I can almost dismiss out of hand. There were some seriously good numerical analysts involved in the writing of those codes, and that is not a mistake they would have made, if it might actually be a "mistake".

Ok, so now I need to return to alternative 2, and discuss why, in fact, a scaling is not necessary, as long as a shift is performed. As long as the lower endpoint is zero, we seriously don't care.

Consider a polynomial that lives on the interval [0,h], and another polynomial that lives on the interval [0,1]. Assuming the two polynomials have exactly the same fundamental shape, except for the scaling, what will they look like? Can we transform one polynomial into the other easily enough?

[0,h] interval:

P_h(x) = a0*x^3 + a1*x^2 + a2*x + a3

[0,1] interval:

P_1(x) = a0*(h*x)^3 + a1*(h*x)^2 + a2*(h*x) + a3

= (a0*h^3)*x^3 + (a1*h^2)*x^2 + (a2*h)*x + a3

So the two polynomials will behave exactly the same! Yes, the coefficients of P_1(x) are different. But look carefully at those coefficients. The interval length merely gets absorbed into the corresponding coefficient from x, which now lives on [0,1]. Therefore we will have no more and no less problems in terms of floating point arithmetic with P_h and P_1.

And I would suggest alternative 2 was the reason, that there is no gain to be had from using a standardized unit interval. You get all of the gain just by making the lower endpoint 0.

[dpb]

@John D'Errico, thanks as always for clear elucidation! It is obvious when actually write it down...as the algebra shows, the amplitude is scaled but the relative change in terms has to be the same. Hadn't recognized the zero shift by itself having that effect before.

[Steven Lord]

Which "builtin tools" are you asking about? Just spline or do you have others in mind? The spline function has been in MATLAB for a long time (it looks like nearly 38 years according to one of the comments):

dbtype 65:70 spline.m

65    %
66    %   See also INTERP1, MAKIMA, PCHIP, PPVAL, MKPP, UNMKPP.
67    
68    %   Carl de Boor 7-2-86
69    %   Copyright 1984-2019 The MathWorks, Inc.
70    

I did a little quick archaeology to determine when the center-and-scale output was added to polyfit etc. It looks like it was in release R12 back in 2000 (about a year before I started at MathWorks), 14 years after spline.m was written by or in conjunction with Prof. de Boor.

[dpb]

I was thinking specifically of both spline and interp1, but also are the options in Curve Fitting TB with fit based on the same algorithms/methods?

[dpb]

RE: your "archeology", I remember polyfit from the early 3.1 release thru R11 before started the self-employed consulting gig (coincidentally also in 2000) and bought the license from former employer and upgraded. Prior to that I had been browbeat/indoctrinated from reading GEP Box that I had built a standardized version around the existing implementation that also included normalizing shape transitions.

Interestingly enough, I hadn't ever taken the time to do the exercise @John D'Errico just illustrated to recognize that the zero shift accomplished the same thing numerically (although doesn't do anything about the shape transformation, of course). I just looked but I do not seem to have an archival copy of that code at hand although I'm sure I could dig it up...