1
$\begingroup$

I am attempting to calculate smoothing spline coefficients based on the description in Reinsch's 1967 paper, but I'm having some trouble. The first derivative is not continuous. Here are the equations I'm attempting to follow:

smoothing spline equations

And this is the MATLAB code. There are some simplifications because all h_sub_i and weight values are equal to 1.

pts = [0      0
  1      1
  2      2
  3      3
  4      2
  5      1
  6      0];

n = length(pts)-1;
p = 10; % closeness of fit

Q = -2*eye(n-1) + diag(ones(1,n-2),1) + diag(ones(1,n-2),-1);
Q(n+1, end) = 0;
Q(n, end) = 1

T = 4/3*eye(n-1) + diag(1/3*ones(1,n-2),1) + diag(1/3*ones(1,n-2),-1)
Ap = transpose(Q)*Q + p*T

c = Ap^-1 * p * transpose(Q) * pts(:, 2);
a = pts(:, 2) - (1/p)*Q*c;
c = [0; c; 0];
d = (c(2:end) - c(1:end-1)) / 3;
b = (a(2:end) - a(1:end-1)) - c(1:end-1) - d;

N = 10000;
breaks = 0:6;
pp = mkpp(breaks, [d(1:n) c(1:n) b(1:n) a(1:n)]);
x = 0:0.01:6;
y = ppval(pp, x);

figure(1)
clf
hold on
grid on
plot(pts(:,1), pts(:,2), '--bo');
plot(x, y, 'r');

The following output is produced:

Q =

    -2     1     0     0     0
     1    -2     1     0     0
     0     1    -2     1     0
     0     0     1    -2     1
     0     0     0     1    -2
     0     0     0     0     1
     0     0     0     0     0


T =

       1.3333      0.33333            0            0            0
      0.33333       1.3333      0.33333            0            0
            0      0.33333       1.3333      0.33333            0
            0            0      0.33333       1.3333      0.33333
            0            0            0      0.33333       1.3333


Ap =

       18.333     -0.66667            1            0            0
     -0.66667       19.333     -0.66667            1            0
            1     -0.66667       19.333     -0.66667            1
            0            1     -0.66667       19.333     -0.66667
            0            0            1     -0.66667       19.333

enter image description here

What is my error?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.