Smoothing Spline Matlab Code (This sparse assembly could be simplified)

SmoothSpline.m

function [mu,m]=SmoothSpline(X,f,p)

% [mu,m]=SmoothSpline(X,f,p)
% 
% This code compute the parameter of "the smoothing natural spline" 
% for the data  (x_i,f_i) with weights {p_i} i.e. the solution of
%
%     \min\Biggl\{ \sum_{i=0}^N p_i\left(\hat{f}(x_i)-f_i \right)^2
%                  +\int_{a}^{b} \bigl|\hat{f}''(x)\bigr|^2dx \Biggr\}
% Input:  
%             X :=  {x_1,...,x_N} 
%             f :=  {f_1,...,f_N} Set of obserbations   
%             p :=  weigths of the method
%
% Output: 
%           mu := {mu_0,...,mu_N} the values of the Smothing Spline at X 
%           m  := Second derivative of the Smothing Spline at X.
% Last modified: March 21, 2018.

N=length(X);  % Number of observations
f=f(:);
p=1./2*p;     %  
P=sparse(1:N,1:N,p,N,N); % P^{-1} Sparse matrix.
h=X(2:end)-X(1:end-1);     % vector with h_{j+1}-h_j.

h2=(1/6)*h(2:end-1);
h3=(1/3)*(X(3:1:end)-X(1:1:end-2)); % H(j) =(1/3)*(h_j+h_{j+1})
CoefR=zeros(3*N-8,1);
CoefR(1:N-2)=1:N-2; CoefR(N-1:2*N-5)=2:N-2; CoefR(2*N-4:end)=1:N-3;
CoefC=zeros(3*N-8,1);
CoefC(1:N-2)=1:N-2; CoefC(N-1:2*N-5)=1:N-3; CoefC(2*N-4:end)=2:N-2;

A=sparse(CoefR,CoefC,[h3 h2 h2],N-2,N-2);
CoefR=repmat((1:N-2)',3,1);
CoefC=zeros(3*N-6,1);
CoefC(1:N-2)=1:N-2; CoefC(N-1:2*N-4)=2:N-1; CoefC(2*N-3:3*N-6)=3:N;
b=[1./h(1:end-1) -(1./h(1:end-1)+1./h(2:end)) 1./h(2:end)];

H=sparse(CoefR,CoefC,b,N-2,N);
h=h(:);
f2=(f(3:end)-f(2:end-1))./(h(2:end))-(f(2:end-1)-f(1:end-2))./h(1:end-1);
f2=f2(:);

m=zeros(N,1); m(2:end-1)=(H*P*H'+A)\f2;
mu=f-P*(H')*m(2:end-1); %  this could be simplified.
end