Skip to content

Instantly share code, notes, and snippets.

@Hugo-Diaz-N
Created October 19, 2018 01:43
Show Gist options
  • Save Hugo-Diaz-N/e2efb1815725eff3e1221872c244e05c to your computer and use it in GitHub Desktop.
Save Hugo-Diaz-N/e2efb1815725eff3e1221872c244e05c to your computer and use it in GitHub Desktop.
Smoothing Spline Matlab Code (This sparse assembly could be simplified)
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment