natural cubic spline interpolation of y-values: how to get derivative of the spline wrt the y-values?

Given a data set with support points x_1,...,x_n and corresponding y-values y_1,...,y_n.

My objective is to create a cubic spline f (with natural boundary conditions) that passes through the y_values. There are, of course, plenty of functions for doing this.

However, for a parameter identification procedure, I have to compute the derivative of the spline f with respect to the y-values -- at arbitrary points within [x1, x_n].

Is there an easy way using built-in functions of Matlab to compute the sensitivities?

[Torsten]

I can understand that you have to compute the derivative of the spline with respect to the parameters, but why with respect to the y-values ?

[SA-W]

Because the values y_1,...,y_n are the parameters of my optimization. My objective is to optimize the values y_1,...,y_n and interpolate between them with cubic splines. Is that clear?

Answer 1

[Bruno Luong]

The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values.

Replace spline command with your function that computes natural spline pp form

x=cumsum(rand(1,10));
y=rand(size(x));



xi=linspace(min(x),max(x),500);
f=ppval(spline(x,y),xi)

f = 1×500

0.8578 0.8215 0.7861 0.7516 0.7180 0.6852 0.6533 0.6223 0.5920 0.5627 0.5341 0.5063 0.4794 0.4532 0.4278 0.4032 0.3794 0.3563 0.3340 0.3124 0.2915 0.2714 0.2519 0.2332 0.2152 0.1978 0.1812 0.1652 0.1498 0.1351

plot(x,y,'or',xi,f,'b') % 

b=eye(length(x));
yd=spline(x,b);
dfdy=ppval(yd,xi); % dfdy(i,j) is the derivative of f(xi(j)) with respect to y(i))

figure
plot(xi,dfdy')

[SA-W]

"The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values. "

I guess you assume f = b_i y_i (sum over i). Then, of course, the derivative of f wrt to y_i is b_i.

However, as for cubic splines, we do not have such a linear combination (correct me if I am wrong) but rather solve a linear system to obtain the coefficents of a local polynomial b_i(x) = a_i +b_i(x-x_i) + c_i(x-x_i)^2 +d_i(x-x_i)^3. The values y_i are hidden in the coefficients a,b,c,d.

[Torsten]

@SA-W

Linearity here means:

spline(y + e_i * h) = spline(y) + h*spline(e_i)

where

e_i = (0,0,...,1,0,0,...,0) (with 1 at the i-th position).

[Bruno Luong]

"Then, of course, the derivative of f wrt to y_i is b_i."

NO

[SA-W]

@Torsten

"Linearity here means:

spline(y + e_i * h) = spline(y) + h*spline(e_i)

Admittedly, I did not think about linearity of this way...

If I do a finite difference calculation, I could directly implement this equation. But how to get an analytical expression for

dSpline/dy(i)

?

[Torsten]

dspline/dy(i) = 
lim (h->0) (spline(y + e_i * h) - spline(y)) / h = 
lim (h->0) (spline(y) + h*spline(e_i)) / h = 
spline(e_i)

Just an Bruno said.

[SA-W]

Does this also apply for the first derivative of the spline?

Say the spline is called f(x) and interpolates the values y1,...,yn at the support points x1,...,xn. Lets define the first derivative s =: df(x)/dx.

Is

ds/dy1 = s(e_i) ?

[Torsten]

The second derivative of a linear function is ... ?

[SA-W]

"The second derivative of a linear function is ... ?"

... of course zero.

Here, f(x) is a cubic spline and I asked for the first derivative s(x) =: df(x)/d(x) which is then a quadratic function. More specifically, I want to evaluate ds/dy1, which is (correct me if I am wrong) also a quadratic function?

[Bruno Luong]

ds/dy1 = d(df/dx)/dy1 = d(df/dy1)/dx.

Use this answer to compute g := df/dy1 as pp form, then use my other answer (with ppder function) to compute dg/dx (which is correctly piecewise quadratic (and NOT quadratic).

[SA-W]

One last thing: Is

spline(y + e_i * h) = spline(y) + h*spline(e_i)

true for all interpolation splines or is it a characteristic of cubic interpolation splines?

[Bruno Luong]

All spline interpolation or least-square fitting with homogeneuous boundary conditions (natural, not a knot, periodic, clamp first derivative to 0) ae linear with respect to y values regardless the degree.

And btw linear mean more generic

spline(y1 + a*y2) = spline(y1) + a*spline(y2)

for vector y1, y2 and scalar a.

[SA-W]

I would like to cite this in a publicaion.

Can you give me a reference where

spline(y1 + a*y2) = spline(y1) + a*spline(y2)

is mentioned in any form?

[Bruno Luong]

You can quote ME.

Seriously this fact is not quite hard to see for people working with spline (and not simply use them).

[SA-W]

Ok. Can you elaborate a little bit more how to "see" this relation? I would really like to understand it and not just using it.

As for me, I know that the coefficients of the piecewise splines are obtained by solving a linear system A*c=b, where the vector b depends (linear?) on the y-values. Was the linear system also your starting point to come up with

spline(y1 + a*y2) = spline(y1) + a*spline(y2)

?

[Bruno Luong]

Or you can use directly the definition of spline:

• Does spline(y1) + a*spline(y2) interpolate y1 + a*y2? Yes.
• Is it a piecewise third order polynomial, Yes.
• Does it have continuous second derivative, Yes
• Does it have zero second derivative at the boundary, Yes.

So spline(y1) + a*spline(y2) is the spline interpolate y1 + a*y2, meaning

spline(y1 + a*y2) = spline(y1) + a*spline(y2)

[SA-W]

That said, if I wanted to calculate

dspline/dy = (spline(y + e_i * h) - spline(y)) / h

by finite differencing,

I could choose an arbitrary value for h as it cancels out anyway, right?

[Bruno Luong]

Right

Answer 2

[Bruno Luong]

In this thread https://fr.mathworks.com/matlabcentral/answers/1894800-how-to-remove-noise-from-curves-and-take-their-derivates?s_tid=srchtitle

you can find my function that compute the derivative of a piecewise polynomiall function (pp), inclusing pp form of the spline functions. This function returns the pp form of the derivative, so you can evaluate using MATLAB ppval.

function ppd = ppder(pp)
ppd = pp;
coefs = ppd.coefs;
n = size(coefs,2);
ppd.coefs = coefs(:,1:n-1).*(n-1:-1:1);
ppd.order = ppd.order-1;
end

[Torsten]

Isn't that computing the derivative with respect to the independent variable x ? OP needs it with respect to the y-values.

[Bruno Luong]

If that is the case:

The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values.

[SA-W]

The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values.

This is not clear to me. If I take the derivative of f (a scalar function) with respect to y_i (a scalar), the result should be a single number and not a vector b_i := (0,0,...,1,0...)...

Also, why is a cubic spline linear to y values? If I change yn to yn+h, f will be more affected at yn than at y1, for instance.

[Bruno Luong]

See my answer below, but

"the result should be a single number"

No the result is a scalar function. If we take f at a given point x then it is a scalar.

"...and not a vector b_i := (0,0,...,1,0...)."

I did not tell the derivative is b, the derivative is the spline interpolating b

"why is a cubic spline linear to y values?"

You clearly missunderstand and male confusion betwen being linear and being a lilnear function

Answer 3

[Torsten]

Why do you want to compute the sensitivities manually ?

Usually, the fitting software computes them using a finite-difference approximation, i.e. by calling your function with

y_i 

and

y_i+h

getting back

f_j(y_1,...,y_i,...,y_n) and f_j(y_1,...,y_i+h,...,y_n)

and approximating

df_j/dy_i = (f_j(y_1,...,y_i+h,...,y_n)-f_j(y_1,...,y_i,...,y_n))/h     

And this would also be my suggestion on how to do it manually if it is really needed.