Direct methods for fitting and evaluating RBFs

This page contains basic matlab code to fit and evaluate an RBF to 3D data using the bi-harmonic basic function, phi(r)=|r|. The direct code is presented here for comparison and clarification of what FarField's FastRBF toolbox does. Essentially, the FastRBF toolbox exhibits the same functionality as the following code except that it is very much faster, numerically more robust, and requires far less memory. These properties become obvious for larger data sets, where the number of data points typically exceeds 2000 (in 3D). Above a few thousand points the solution and evaluation of an RBF by direct methods becomes impossible on most machines. 
% Demonstration file for a simple RBF biharmonic spline 
% fit in 3D
%
%
n=100;
cent=rand(100,3);
f=rand(100,1);
coeff=fitit(cent,f);
fprintf('Value to be interpolated  %15.10e \n',f(1));
fprintf('Value of RBF interpolant  %15.10e \n',...
       eval_direct(cent,coeff,cent(1,:)));

fitit.m - Direct fit of an RBF

Fitting an RBF invloves setting up and solving the following linear system. 
function [coeff]=fitit(cent,f)
% Function to find the radial basis function s consisting
% of a linear plus a sum of shifts of \phi( | x | )
% interpolating to data
% f_i at the point cent(i,:) for 1 \leq i \leq n.
%
% Syntax [coeff]=fitit(cent,f)
%
% Input
%
%    cent    n by dim array of centres
%    f       n by 1 vector of values at centres
%
% Output
%    coeff   (n +3) by 1 vector of coefficients with the 
%            coefficients of 1, x and y (2D) or 1, x, y 
%            and z (3D) last.
%
A=direct(cent);
%
[n dim] =size(cent);
f=f(:);                % Make sure its a column
f=[f ;zeros(dim+1,1)]; % add zeros for the polynomial part
                       % at the end of the column
coeff=A\f;

direct.m - Forming the interpolation matrix

function [A]=direct(cent)
% Form the (n+dim+1)*(n+dim+1) matrix corresponding to 
% RBF interpolation with basic function phi and linears.
%
% The centres are assumed given in the n by dim array cent.
% phi is assumed given as a function of r. It is coded in 
% the Matlab function phi. m
%
% Syntax  [A]=direct(cent)
%
% Input          
%         cent  n*dim array  coordinates of centers
%               of reals     
% Output  A     (n+dim+1)*   Symmetric matrix of the linear
%               (n+dim+1)    system obtained if we solve
%               array of     for the radial basis function
%                            interpolant directly.
%
% Write the matrix A in the form
%             B    P
%     A   =  
%             P^t  O
% where P is the polynomial bit.
%
[n dim]=size(cent);
A=zeros(n,n);
for i=1:n
    for j=1:i
        r=norm(cent(i,:)-cent(j,:));
        temp=phi(r);
        A(i,j)=temp;
        A(j,i)=temp;
    end
end
%
% Now the polynomial part
%
P=[ones(n,1) cent];
A = [ A      P
      P' zeros(dim+1,dim+1)];

phi.m - Basic function

Here we consider the bi-harmonic basic function in 3D, phi(r)=|r|. 
function u=phi(r)
%
% Function phi of r
%
% Syntax [u] = phi(r)
%
% Remember if using something like the thinplate spline
 
% in 2D you will need to test for r nonzero before 
% taking the log.
u =r;

eval_direct.m - Direct evaluation

Evaluating the RBF $s(\vec{x})$. 
function s=eval_direct(cent,coeff,u)
%
% Function which directly evaluates the RBF spline
%
% s(u)=linear poly +\sum_i=1^n coeff(i)*phi(u-cent(i,:) 
%
% at the point u where the linear polynomial is given  
% in terms of its monomial cofficients. 
%
% Syntax     s=eval_direct(cent,coeff,u)
%
% Inputs
%    cent   n by dim array of coordinates for the centres
%    coeff  n+dim +1 vector of coefficients with the linear
%           polynomial part last.
%    u      row vector    point at which to evaluate
%
% Output
%    s    Value of the RBF at position u
%
% Code is written for clarity rather than Matlab efficiency
[n dim] = size(cent);
s=0;
for j=1:n
   s = s + coeff(j)*phi(  norm(u - cent(j,:) )   );
end
% Now the linear polynomial bit
s=s+coeff(n+1);               % The constant
for i=1:dim
   s=s+coeff(i+n+1)*u(i);     % The various components of u
end

FastRBF FAQ