express basis spline derivative in terms of interpolation values

x = [1 2 3 4 5];
y = randn(size(x));
f = spapi(5,x,y);
fdd = fnder(f, 2);

If we evaluate the second derivative at a given point x* as

xeval = 3.5;
val = fnval(fdd, xeval)

is there a way to express "val" as a linear combination of the interpolation values as

val = y(1)*c1 + ... + y(5)*c5

Then, the task is to find the coefficients ci.

Can this be done with symbolic differentiation or other techniques? Or is there an easy analytical representation for that?

Thank you!

Answer 1

[Bruno Luong]

x = [1 2 3 4 5];
y = randn(size(x));
k = 5;

f = spapi(k,x,y);
B = spapi(k,x,eye(length(x)));

fdd = fnder(f, 2);
Bdd = fnder(B,2);

xeval = 3.5;
val = fnval(fdd, xeval)

val = 1.5349

c = fnval(Bdd, xeval);
vcomb = y * c

vcomb = 1.5349

err = val-vcomb

err = 8.8818e-16

[SA-W]

Thank you! The explanation here is that the spline is linear w.r.t y?

[Bruno Luong]

Well we discussed many time. I repeat myself

• spline function s is linear wrt to y; remember how to compute the Basis of s wrf to y using eye(...)?
• s'' is linear wrt s

Therefore

• s'' is linear wrt y (thus your question). Just chaining the two above you'll get the basis of s"" wrt to y. It seems not difficult to make this deduction.

[SA-W]

Makes sense.

May I ask a follow up-question for the case spline order k=5, hence second derivative is piecewise quadratic.

When I choose

% x: points
% y: values
xeval = x;
A = fnval(Bdd, xeval)';
A*y >= 0

the second derivative will be postive at the breaks but can be negative between the breaks. To approximately have it postive on the entire domain, I could choose xeval = linspace(min(x), max(x), n) with n very large. Is there a better solution how to choose xeval sucht that f''>= 0 on the entire interpolation domain? I do not want to create too many constraints (number of rows of A =: n), hence keep n as small as possible but as large as needed to encode f''>=0

[Bruno Luong]

I wouldn' formalize the constraint s''(x) >= 0 for all x in term of not sufficient) condition s''(ti) >= 0 where ti is the knots, as you have rightly noticed.

Rather in the B-spline form if I express:

s''(x) = sum_i gamma_i N_{i,2}

where N_{i,2} is B-spline basis of order 3 (quadratic spline) and gamma_i is coeffiecients of the second derivative s''.

We know N_{i,2}(x) (as the case of all B-spline basis) is positive and sum to unity, the sufficient condition for s''(x) >= 0 is to simply require gamme_i >= 0 for all i.

The relationshif of y -> gamma is linear as we have noted about the linearity of spline wrt to y. A little of spline command and matrix algebra you can the express gamma = G*y, where G is some matrix.

That what I have implemented in my BSFK package, and so far it seems working fine and reliable as the few examples I gave you in this early thread.

[Bruno Luong]

I'll repost the example here with fifth order spline

x = [3, 3.8125, 4.6250, 5.4375, 6.25, 7.0625, 7.8750, 8.6875, 9.5, 10.3125, 16];
y = [0, 0.0198463, 0.0397084, 0.0596019, 0.0795445, 0.0995199, 0.119908, 0.140298, 0.160773, 0.181772, 1.02478];
% FEX https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation?s_tid=srchtitle_BSFK_1
opt = struct('KnotRemoval','none','sigma',1e-10,'shape',struct('p',2,'lo',0,'up',inf));
k = 5; % spline order
pp = BSFK(x,y,k,20,[],opt);
% prepare graphical data
xi = linspace(min(x),max(x),1025);
yi = ppval(pp,xi);
xb = pp.breaks(2);
yb = ppval(pp,xb);
% Check if approximation is close to interpolation
norm(y-ppval(pp,x),'inf')/norm(y,'inf') % 2.6938e-05
% Evaluate secpn derivative
fdd = ppder(ppder(pp));
yddi = ppval(fdd,xi);
all(yddi >= 0) % true
figure(1)
clf
plot(x, y, 'or')
hold on
plot(xi, yi, 'b')
xlabel('x')
ylabel('y')
yyaxis right
plot(xi, yddi) % second derivative is postive
ylabel('y''''')
legend('data', 'spline interpolation', 'second derivative')

[SA-W]

Thats indeed way simpler because gamme_i >= 0 for all i are bound constraints and no inequality constraints.

But it changes my method a lot and for now, I would prefer a solution to express s''(x)>=0 in terms of y, hence in pp-form. Your answer shows how to impose this constraint at discrete grid points xi, but we have no control what happens between the xi's. Simply imposing the constraint at many points will increase the condition number of A and may case numerical troubles.

I just tried

xeval = linspace(min(x), max(x), n) % n >> numel(x)
A = fnval(Bdd, xeval)';
A = licols(A')'

which seems to work but I am not sure how this performs generally.

What do you think of this solution?

And can you think of alternatives?

[Bruno Luong]

Honestly I don't like your method because this do not ensure s''(x) >= 0 with certainty as you notices that make your matrix ill-condition when you increase the density and use li to remove the points (column of A" or row of A). IMO You better keep the matrix without reduce with licols. It is not clear to me the role of repeat knots on the boundary in your method

"And can you think of alternatives?"

Why keep asking my opinion if you want to stick with your method?

[Bruno Luong]

BTW if you use cubic spline, the second derivative is piecewise linear and impose positivity s''(ti)>=0 at the knots ti (better than xi) will ensure s''(x) >= 0 equavilently.

In short the usage of fifth order spline makes your constraints more difficult to deal with.

[SA-W]

In short the usage of fifth order spline makes your constraints difficult to deal with.

I know, but a piecewise linear second derivative is not "smooth" enough for my application. Thats what I learned so far.

The relationshif of y -> gamma is linear as we have noted about the linearity of spline wrt to y. A little of spline command and matrix algebra you can the express gamma = G*y, where G is some matrix.

So in this notation, G is the matrix of linear inequality constraints (G*y >= 0).

Can you help me to construct G?

[SA-W]

Here is some code aiming to compute G:

k = 5;
B = spapi(k,x,eye(numel(x)));
Bdd = fnder(B, 2);
A = fnval(Bdd, x)';

t2 = aptknt(x, k-2); % knots for spline s''(x) = \sum_i N_{i,2}(x) gamma_i
C = spcol(t2, k-2, x); % C(i,j) = N_{i,2} evaluated @ x(j)

G = inv(C)*A;
% G*y >= 0 --> s''(x)>=0 on interpolation domain

Let s''(x) = sum_i gamma_i N_{i,2} be the second derivative (piecewise quadratic). In your answer, you mentioned the equation f''(x) = Bdd * y. Obviously, they must match: s''(x) = f''(x) \forall x. Theres are numel(x) unknows gamma_i, hence we can s''(xj) = f''(xj) --> N_{i,2}(xj) gamma_i = Bdd(xj) * y_i for i,j=1,...,numel(x).

Is that correct?

[Bruno Luong]

Can you help me to construct G?

It should be like this

x = [1 2 3 4 5];

k = 5;
B = spapi(k,x,eye(length(x)));
Bdd = fnder(B,2);
G = Bdd.coefs'

G = 3x5

2.9167 -8.6667 9.5000 -4.6667 0.9167 -2.0833 9.3333 -14.5000 9.3333 -2.0833 0.9167 -4.6667 9.5000 -8.6667 2.9167

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

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

[SA-W]

Nice, thanks!

Maybe you can comment on my approach below:

k = 5;
x = [1 2 3 4 5];
B = spapi(k, x, eye(numel(x)));
Bdd = fnder(B, 2);
k2 = 3;
t2 = aptknt(x,k2);
C = spcol(t2,k2,x);
A = fnval(Bdd, x)';
G = inv(C)*A;
rank(G); % 3

If I am not mistaken, you make the ansatz s''(x) \sum_i N_{i,2}(x) gamma_i where i=1,2,3 resulting in a 3x5 matrix. I make the same ansatz but with i=1,2,3,4,5 resulting in a 5x5 matrix with rank(G) = 3 and cond(G) close to infinity. These numbers tell me I have too many constraints (rows) and your 3x5 matrix is more appropriate. That said, my matrix is also correct to enforce s''(x) = G*y >= 0 but contains redundant information, right?

[Bruno Luong]

I don't look carefully what you did. gamma(:) is 3 x 1 so B should be 3 x 5 matrix.

Why make those spagetti with x?

[SA-W]

Makes sense.

If I take your code with

y = [0 0.00634 0.0148 0.0246 0.0356];

we obtain G*y =

G*y =

0.003486666666667
0.000006666666667
0.001646666666667

The first and last value is the second derivative at the left and right point, respectively. However, 0.000006666666667 does not match the second derivative at some point. So can we interpret the meaning of the, in this case, second row of G?

[Bruno Luong]

Gamma = B*y is interpreted the coefficients of the second derivative function in B-spline basis, not the derivative ifself. The coefficients are people usually call in spline context control point (here 1D)

[SA-W]

Oh yes, I mixed some terms.

Last point:

k = 5;
x = [1 2 3 4 5];
y = [0 0.00634 0.0148 0.0246 0.0356];
f = spapi(k,x,y);
fdd = fnder(f, 2);
fextrap = fn2fm(fnxtr(f, k), 'B-');
fddextrap = fn2fm(fnxtr(fdd, k), 'B-')

fddextrap = struct with fields:
      form: 'B-'
     knots: [0 0 0 6 6 6]
     coefs: [0.0055 -0.0016 0.0028]
    number: 3
     order: 3
       dim: 1

Suppose I want to extrapolate the spline using the same order as shown above. Here, we see that gamma coefficients can be negative although the second derivative is positive eveywhere. Hence, gamma_i >= 0 no longer holds.

Did I make a mistake or is not so straightforward to account to express s''>=0 for the extrapolating spline?

[SA-W]

k = 5;
x = [3 4 5 6 7 8 9 10 11 12 16];
y = [0.0000000000000
0.0000578540425
0.0004881270575
0.0009749918603
0.0015958783882
0.0023269066095
0.0031921431120
0.0041236074503
0.0051080185669
0.0071627232461
0.2156297589267];

f = spapi(k,x,y);
fdd = fnder(f, 2)

fdd = struct with fields:
      form: 'B-'
     knots: [3 3 3 5.5000 6.5000 7.5000 8.5000 9.5000 10.5000 16 16 16]
     coefs: [0.0012 -2.9217e-04 2.3087e-04 4.0579e-05 2.0349e-04 1.7634e-05 7.8003e-05 -3.6568e-04 0.0823]
    number: 9
     order: 3
       dim: 1

In this example (without extrapolation), the second derivative is positive everywhere but some of the gamma coefficients are negative. Why is that?

[Bruno Luong]

Why not?

Read very carefully what I wrote

"We know N_{i,2}(x) (as the case of all B-spline basis) is positive and sum to unity, the sufficient condition for s''(x) >= 0 is to simply require gamma_i >= 0 for all i."

There is nothing in your last two questions contracdict with what I wrote.

Try to backup your expectation in the questions: why coefs should be >= 0, did I wrote something that makes you think it should?

Sigh, you always make false deduction.

[SA-W]

I am afraid if I misunderstood you again, sorry!

But my understanding is that

gamma = fdd.coefs

and from

gamma_i >= 0 follows s''(x) >= 0 (sufficient condition).

Do I deduce this correct?

In the example I posted, we have the case s''(x)>= 0 but gamma_i < 0. Does it mean s''(x)>=0 is possible although some gamma_i < 0 ?

[Bruno Luong]

"gamme_i >= 0" is the sufficient condition for s" >= 0 meaning

gamma_i >= 0 implies s''(x) >= 0, or

if gamma_i >= 0 then s''(x) >= 0

Not the opposite that I NEVER state.

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

[SA-W]

Yes, it is clear now.

k = 5;
x = [1 2 3 4 5 6 7 8];
y = randn(size(x));
f = spapi(k,x,y);
fdd = fnder(f, 2);

k = f.order;
t = f.knots;
fi = spmak(t(1:k+1), 1);
fnplt(fi);
hold on;
for i=2:numel(f.coefs)
   fi = spmak(t(i:i+k), 1);
   fnplt(fi);
end

figure(2);
k = fdd.order;
t = fdd.knots;
fi = spmak(t(1:k+1), 1);
fnplt(fi);
hold on;
for i=2:numel(fdd.coefs)
   fi = spmak(t(i:i+k), 1);
   fnplt(fi);
end

Is this code correct to plot the basis functions N{i,5} (s(x) = \sum_i N{i,5} p_i) and N{i,3} (s''(x) = \sum_i N{i,3} gamma_i)

[Bruno Luong]

Another way of plotting basis B-spline, no need to create dummy y, spapi, and fnder.

Also note that my "fi" are defined on the whole interval covering x, not just on their respective supports [t(i) t(i+k)]

k = 5;
x = [1 2 3 4 5 6 7 8];
t = aptknt(x,k);
n = length(t) - k;
figure
hold on
for i = 1:n
    fnplt(spmak(t, accumarray(i, 1, [n 1])'));
end
% Check if they sum to 1
fnplt(spmak(t, ones(1,n)),'k-.');

d = 2; % derivative order
k = k-d;
t = t(1+d:end-d);
n = length(t) - k;
figure
hold on
for i = 1:n
    fnplt(spmak(t, accumarray(i, 1, [n 1])'));
end
% Check if they sum to 1
fnplt(spmak(t, ones(1,n)),'k-.');

[SA-W]

Thanks! I appreciate your expertise and support, but feel free to stop the conversation.

I was just wondering about the construction of the second derivative basis splines. Given s(x) = \sum_i N_{i,k}(x) P_i with i=1,...,n, it is easy to understand that the basis functions N_{i,k}>=0. With spline order k, the knot sequence has n+k elements. So far so good.

Given that control points P_i are constant, we can compute the second derivative by differenting the basis functions twice: s''(x) = \sum_i N''_{i,k}(x) P_i with i=1,...,n, resulting in a spline of order k-2, hence, n+k-2 elements of the knot vector. Obviously, N''_{i,k}(x) >=0 is false.

Your code snippet and fnder(f,2) produce s''(x) = \sum_i N_{i,k-2} gamma_i, where i=1,...,n-2. Here, N_{i,k-2}>=0 . I do not see how this version of the second derivative is related to the orignal spline s(x).

Is my version (s''(x) = \sum_i N''_{i,k}(x) P_i) just a different way of expressing the second derivative or is it completely wrong?

[Bruno Luong]

" I do not see how this version of the second derivative is related to the orignal spline s(x)."

N''_{i,k}(x) is NOT a basis since it is not independent, not positive and and it is not B-spline basis standard as defined by De Boor and friends.

To express in B-spline basis you need to read about the spline derivative literature, for example here

https://public.vrac.iastate.edu/~oliver/courses/me625/week5b.pdf

Such relation is surely implemented under the hood in MATLAB fnder function.

[SA-W]

Such relation is surely implemented under the hood in MATLAB fnder function.

Indeed. I double-checked that coefficients of first derivative are computed according to the equation given on page 203 (Qi = ...) in the linked document.

But I think my approach was not wrong: In the document, they start with s(x) = \sum_i N_{i,k}(x) P_i, then differentiating s'(x) = \sum_i N'_{i,k}(x) P_i, and do some manipulations do arrive at s'(x) = \sum_i N_{i,k-1}(x) Q_i.

[Bruno Luong]

"But I think my approach was not wrong"

No it's not. But you seem to question about why s'' can be expressed as apparently by two different ways, and how to transform the first one to the second one.

I hope now it is clear for you.

[SA-W]

To conclude this here:

B = spapi(k, x, eye(numel(x)));
Bdd = fnder(B, 2);
G = Bdd.coefs';
Bextrap = fn2fm(fnxtr(B, k), 'B-');
Bddextrap = fn2fm(fnxtr(Bdd, k), 'B-');
Gextrap = Bddextrap.coefs';

When I want to extrapolate the spline, I can use Gextrap to have s''(x)>=0 outside the interpolation domain.

[Bruno Luong]

Note that you can compute fddextrap differently, but be aware if the extrapolation order changes, your method returns the extrapolation low order spline that is not the second derivative of the extrapolation.

k = 5;
x = [1 2 3 4 5 6 7 8];
y = randn(size(x));
f = spapi(k, x, y);
fdd = fnder(f, 2);
fextrap = fn2fm(fnxtr(f, k), 'B-');

% These three are equal
fddextrap1 = fn2fm(fnxtr(fdd, k), 'B-')

fddextrap1 = struct with fields:
      form: 'B-'
     knots: [0 0 0 3.5000 4.5000 5.5000 9 9 9]
     coefs: [-26.0383 10.9577 -3.7168 0.6648 2.6367 -10.3078]
    number: 6
     order: 3
       dim: 1

fddextrap2 = fn2fm(fnxtr(fdd, k-2), 'B-')

fddextrap2 = struct with fields:
      form: 'B-'
     knots: [0 0 0 3.5000 4.5000 5.5000 9 9 9]
     coefs: [-26.0383 10.9577 -3.7168 0.6648 2.6367 -10.3078]
    number: 6
     order: 3
       dim: 1

fddextrap3 = fnder(fextrap,2)

fddextrap3 = struct with fields:
      form: 'B-'
     knots: [0 0 0 3.5000 4.5000 5.5000 9 9 9]
     coefs: [-26.0383 10.9577 -3.7168 0.6648 2.6367 -10.3078]
    number: 6
     order: 3
       dim: 1

fextrap_k2 = fn2fm(fnxtr(f, 2), 'B-');

% These however are NOT equal
fddextrap_k2_1 = fn2fm(fnxtr(fdd, 2), 'B-')

fddextrap_k2_1 = struct with fields:
      form: 'B-'
     knots: [0 0 0 1 3.5000 4.5000 5.5000 8 9 9 9]
     coefs: [-22.0865 -15.4680 7.6967 -3.7168 0.6648 2.1985 -6.6094 -9.1259]
    number: 8
     order: 3
       dim: 1

fddextrap_k2_2 = fn2fm(fnxtr(fdd, 0), 'B-')

fddextrap_k2_2 = struct with fields:
      form: 'B-'
     knots: [0 0 0 1 1 1 3.5000 4.5000 5.5000 8 8 8 9 9 9]
     coefs: [0 0 0 -8.8495 7.6967 -3.7168 0.6648 2.1985 -4.0928 0 0 0]
    number: 12
     order: 3
       dim: 1

fddextrap2_k2_3 = fnder(fextrap_k2,2)

fddextrap2_k2_3 = struct with fields:
      form: 'B-'
     knots: [0 0 0 1 1 1 3.5000 4.5000 5.5000 8 8 8 9 9 9]
     coefs: [0 0 0 -8.8495 7.6967 -3.7168 0.6648 2.1985 -4.0928 1.3323e-15 0 -1.3323e-15]
    number: 12
     order: 3
       dim: 1

[SA-W]

Thanks for pointing out!

B = spapi(k, x, eye(numel(x)));
Bdd = fnder(B, 2);
G = Bdd.coefs';
Bextrap = fn2fm(fnxtr(B, k), 'B-');
Bddextrap = fn2fm(fnxtr(Bdd, k), 'B-');
Gextrap = Bddextrap.coefs';

I only consider the case that the extrapolation order is the same as interpolation order (continue first and last spline if we think in pp-form)

So given that, Gextrap as calculated is correct to implement Gextrap * y(:) >= 0, is it?

[Bruno Luong]

"So given that, Gextrap as calculated is correct to implement Gextrap * y(:) >= 0, is it?"

Yes. If the goal is to interpoate the (unknown) data x/y such that s(x) is convex on larger domain.

But keep in mind that extrapolation, especially with original/high order is an unstable operation, so dangeruous. It depends on what you do after, it might or might not matter.

The picture in the doc file https://www.mathworks.com/help/curvefit/fnxtr.html

illustrate the problem of extrapolate when using the original spline beyond the domain (blue curve). Lower order to 2 will make extrapolation more stable (red).

BTW I did not fully fully understand this picture since the spline interpolate x -> x.^3 in 21 points and I can't see the reason why the extrapolation (blue curve) represents different function (curved concavely) outside the domain.

[SA-W]

BTW I did not fully fully understand this picture since the spline interpolate x -> x.^3 in 21 points and I can't see the reason why the extrapolation (blue curve) represents different function (curved concavely) outside the domain.

I think its because the extrapolation spline is order 2 (linear) while the original spline is order 4 (cubic).

Yes. If the goal is to interpoate the (unknown) data x/y such that s(x) is convex on larger domain.

I am aware of that, but I have to offer it for a small percentage of evaluations.

But I just think whether I need the extrapolation spline. x(1) is an lower bound for the x variable, hence I never need to extrapolate on the left end but only on the right end x(end).

To make s(x) convex beyond x(end) is is not simply sufficient to implement s''''(x(end))>0 (fourth derivative). That is, we have Bdd.coefs' to make it convex on basic interval + one more line fnval(Bdddd, x(end))' to have positive fourth derivative at the end. In code:

k = 5
x = [1 2 3 4 5];
B = spapi(k, x, eye(numel(x)));
Bdd = fnder(B, 2);
G = Bdd.coefs';
% add s''''(x)>=0
Bdddd = fnder(B, 4);
G = [G; fnval(Bdddd, x(end))'];

Does that make sense or would you not do it that way?

[Bruno Luong]

I think its because the extrapolation spline is order 2 (linear) while the original spline is order 4 (cubic).

Not that but it seems csaps does non intuitive thing to "smooth" the spline (minimize the second derivative L^2 norm) and resulting strange extrapolation instead of extrapolate with x^3 as with straight forward spline() command:

x = linspace(0,1,21);
y = x.^3;
k = 4;
f = csaps(x,y); 
extrSpline = fnxtr(f,k);


s = spline(x,y); 
extrSpline = fnxtr(s,k);


figure
for d=0:3
    fd = fnder(f, d);
    subplot(2,2,d+1);
    fnplt(fd,[-.5 1.4])
    sd = fnder(s, d);
    hold on
    fnplt(sd,[-.5 1.4])
    legend('csaps','spline','location','best')
    xline(max(x), 'DisplayName', 'right bound')
    title(sprintf('%d order derivative', d))
end

Does that make sense or would you not do it that way?

Honestly I dont understand what you want to achieve your code. Why do you put positive constraint on fdddd(x(end)) and not fddd(x(end))? I don't know what is your goal.

But never mind I don't really want to discuss those details, you seem to know what you are doing.

[SA-W]

Honestly I dont understand what you want to achieve your code. Why do you put positive constraint on fdddd(x(end)) and not fddd(x(end))? I don't know what is your goal

Oh, you are right. To have s''(x)>0 for x>x(end), I thought the fourth derivative must be positive at x(end). But s'''(x(end)) is sufficient, right?

(I will not ask further after this and say thanks a lot for helping!! :-))

[Bruno Luong]

You might need to impose both conditions of third and fourth derivative on x(end).