Spline interpolation of multivariate data

Hi,

I need some guidance on how to perform cubic spline interpolation of data. The process I'm trying to model has three inputs (independent variables) and three outputs (dependent variables). As I do not know much about splines, I need some clarification on following questions:

1) Among tensor product splines, is there a method that maps R^3 -> R^3 (similar to regression models), or each dependent variable requires a separate model?

2) As examples I found all use regular grids, how can cubic spline interpolation of non uniformly spaced points in R^3 be implemented in MATLAB?

Thanks, Davor

[Shrinivas Iyengar]

Hi @Davor, did you ever manage to find an answer to this?

[Davor]

Hi @Shrinivas Iyengar, I used radial basis functions. Here's my implementation. I won't explain this in detail now since I'm home on vacation with 2 week old newborn. Let me know if you have questions and I'll try to find time to provide explanations.

classdef rbf
    %RBF Radial basis function
    %   Detailed explanation goes here
    
    properties
        lambda
        centers
        kernel
        epsilon
        Xmin
        Xmax
        Ymin
        Ymax
    end
    
    methods
        function obj=rbf(in, out, kernel, epsilon)
            X=in;
            Y=out;
            obj.kernel=kernel;
            obj.epsilon=epsilon;
            obj.Xmin=min(min(X));
            Xs=X-obj.Xmin;
            obj.Xmax=max(max(Xs));
            Xs=Xs/obj.Xmax;
            obj.centers=Xs;
            obj.Ymin=min(min(Y));
            Ys=Y-obj.Ymin;
            obj.Ymax=max(max(Ys));
            Ys=Ys/obj.Ymax;
          %Vectorize code. Set input points as centers. As centers, they
          %are in rows, repaeted in columns.
            centersmat=repmat(obj.centers,size(X,1),1);
            centersmat=reshape(centersmat,size(obj.centers,1),size(X,1),size(X,2));
            %Setting input points (now as samples, not centers!) as columns
            %requires loop. Setting centers as rows and samples as columns
            %enables subtraction of each point (sample) from each point (center).
            for i=1:1:size(X,2)
                Xsmat(1:size(obj.centers,1),1:size(X,1),i)=repmat(Xs(:,i)',size(obj.centers,1),1);
            end
            A=centersmat-Xsmat;
            A=A.^2;
            radiusessquare=sum(A,3);
            switch kernel
                case 'gaussian'
                    %Gaussian
                    A=exp(-radiusessquare*obj.epsilon);
            %Thin plate spline r^2*log(r)
                %A=radiusessquare.*log(sqrt(radiusessquare));
                %A(isnan(A))=0;
                case 'polyharmonic'
            %Polyharmonic spline r^4*log(r)
                %A=radiusessquare.^2.*log(sqrt(radiusessquare)*obj.epsilon);
                %Proba potencija 3
                A=radiusessquare.*sqrt(radiusessquare);
                A(isnan(A))=0;
            %Polyharmonic spline r^6*log(r)
                %A=radiusessquare.^3.*log(sqrt(radiusessquare));
                %A(isnan(A))=0;
            end
            %obj.lambda=pinv(A)*Ys;
            obj.lambda=A\Ys;
        end
        function out = eval(obj, in)
            X=in;
            Xs=X-obj.Xmin;
            Xs=Xs/obj.Xmax;
            %"out" is calculated, predicted by model. Initialize "out".
            out=zeros(size(Xs,1),size(obj.lambda,2));

        %Vectorize code. Centers are set as rows and repeated in columns.
        %Input points are set as columns and repeated in rows. This enables
        %subtraction of each point from each center For large number of
        %inputs, this leads to memory issues. Thus, for many inputs,
        %meshgrid is limited to 100000000 entries and inputs broken to parts.

            if size(X,1)*size(obj.centers,1)>1000000;%20000000
                no_cols=round(20000000/size(obj.centers,1));
                for j=1:no_cols:size(X,1)
                    %Setting input points as columns requires loop.
                    for i=1:1:size(X,2)
                        if j+no_cols-1<=size(X,1)
                            Xsmat(1:size(obj.centers,1),1:no_cols,i)=repmat(Xs(j:j+no_cols-1,i)',size(obj.centers,1),1);
                        else
                            Xsmat(1:size(obj.centers,1),1:(size(X,1)-j+1),i)=repmat(Xs(j:size(X,1),i)',size(obj.centers,1),1);
                        end
                    end
                    %Calculate part of solutions and append to output
                    centersmat=repmat(obj.centers,size(Xsmat,2),1);
                    centersmat=reshape(centersmat,size(obj.centers,1),size(Xsmat,2),size(X,2));
                    diff=centersmat-Xsmat;
                    diff=diff.^2;
                    radiusessquare=sum(diff,3);
                    switch obj.kernel
                        case 'gaussian'
                            %Gaussian
                            phi=exp(-radiusessquare*obj.epsilon);
                    %Thin plate spline r^2*log(r)
                        %phi=radiusessquare.*log(sqrt(radiusessquare));
                        %phi(isnan(phi))=0;
                        case 'polyharmonic'
                            %Polyharmonic spline r^4*log(r)
                            %phi=radiusessquare.^2.*log(sqrt(radiusessquare)*obj.epsilon);
                            %Proba potencija 3
                            phi=radiusessquare.*sqrt(radiusessquare);
                            phi(isnan(phi))=0;
                    %Polyharmonic spline r^6*log(r)
                        %phi=radiusessquare.^3.*log(sqrt(radiusessquare));
                        %phi(isnan(phi))=0;
                    end
                    phi=phi';
                    for i=1:1:size(obj.lambda,2)
                         out(j:j+size(Xsmat,2)-1,i)=sum(bsxfun(@times,phi,obj.lambda(:,i)'),2); 
                    end
                    clear Xsmat;
                end
                out=out*obj.Ymax;
                out=out+obj.Ymin;
            else
                centersmat=repmat(obj.centers,size(X,1),1);
                centersmat=reshape(centersmat,size(obj.centers,1),size(X,1),size(X,2));
                %Setting input points as columns requires loop.
                for i=1:1:size(X,2)
                    Xsmat(1:size(obj.centers,1),1:size(X,1),i)=repmat(Xs(:,i)',size(obj.centers,1),1);
                end
                diff=centersmat-Xsmat;
                diff=diff.^2;
                radiusessquare=sum(diff,3);
                switch obj.kernel
                    case 'gaussian'
                    %Gaussian
                    phi=exp(-radiusessquare*obj.epsilon);
                %Thin plate spline r^2*log(r)
                    %phi=radiusessquare.*log(sqrt(radiusessquare));
                    %phi(isnan(phi))=0;
                    case 'polyharmonic'
                %Polyharmonic spline r^4*log(r)
                        %phi=radiusessquare.^2.*log(sqrt(radiusessquare)*obj.epsilon);
                        %Proba potencija 3
                        phi=radiusessquare.*sqrt(radiusessquare);
                        phi(isnan(phi))=0;
                %Polyharmonic spline r^6*log(r)
                        %phi=radiusessquare.^3.*log(sqrt(radiusessquare));
                        %phi(isnan(phi))=0;
                end
                phi=phi';
                for i=1:1:size(obj.lambda,2)
                   out(:,i)=sum(bsxfun(@times,phi,obj.lambda(:,i)'),2); 
                end
                out=out*obj.Ymax;
                out=out+obj.Ymin;
            end
        end
    end
    
end