Interpolation using compactly supported radial basis function

Question

I have been struggling for two days with the following problem.

I would like to do a $d$-dimensional interpolation over some data. I tried first to use polyharmonic splines, but when the size of data increase, the inversion of matrices is problematic because they are almost singular (large condition number). Therefore, I tried to use compactly supported radial basis function (CSRBF).

What I cannot understand is: why using exactly the same algorithm (see below), it works with a Thin-Plate Spline function (or any other polyharmonic basis function) but it fails completely with any CSRBF.

Here I provide a minimal partially working example that implements two working functions (see wiki)

• $r^2\ln(r)$
• $r^3$

and two non-working functions (see page 39 of the pdf)

• $(1-r)^2$
• $(1-r)^4_+(4r+1)$

here is the code, I tried to simplify it as much as I could, so it is not efficient but you can play with it and see my comments :

#include <vector>
#include <iostream>
#include <random>

unsigned int N(10);
std::vector<double> x;
std::vector<double> y;
std::vector<double> z;
std::vector<double> weights;

extern "C" void dsytrf_(char const& uplo, unsigned int const& n, double *m, unsigned int const& lda, int *ipiv, double *work, int const& lwork, int& info);
extern "C" void dsytri_(char const& uplo, unsigned int const& n, double *m, unsigned int const& lda, int const *ipiv, double *work, int& info);

class Matrix{
    public:
        Matrix(unsigned int L):L_(L),mat_(new double[L_*L_]){}
        ~Matrix(){ delete[]  mat_; }
        double const& operator()(unsigned int const& i, unsigned int const& j) const { return mat_[i+j*L_]; }
        double& operator()(unsigned int const& i, unsigned int const& j) { return mat_[i+j*L_]; }
        /*method that calls lapack routine to invert a symmetric matrix (only
         *need the upper triangular part)*/
        void inv(){
            int info(1);
            int* ipiv(new int[L_]);
            double* work(new double[L_]);
            dsytrf_('U',L_,mat_,L_, ipiv, work,L_,info);
            delete[] work;
            work = new double[L_];
            dsytri_('U',L_,mat_,L_,ipiv,work,info);
            delete[] work;
            delete[] ipiv;
            for(unsigned int i(0);i<L_;i++)
                for(unsigned int j(i+1);j<L_;j++)
                    (*this)(j,i) = (*this)(i,j);
        }

    private:
        unsigned int L_;
        double* mat_;
};
/*everything above this line is irrelevent to my question*/

/*with the first or second f, I have the expected behaviour, but if I use the
 *third or fourth definition of f, it fails*/
double f(double r){ return r*r*r; }
//double f(double r){ return r*(r<1?log(pow(r,r)):r*log(r)); }
//double f(double r){ return (r<1.0?(1-r)*(1-r):0.0); }
//double f(double r){ return (r<1.0?(1-r)*(1-r)*(1-r)*(1-r)*(4*r+1):0.0); }

/*once the weights have been computed, this function returns the interpolated
 * z-value. if (a,b) is one of the known points (x,y), is should exactly the
 * same value z*/
double extrapolate(double a, double b){
    double tmp(0);
    for(unsigned int i(0);i<N;i++)
        tmp += weights[i]*f(sqrt( (a-x[i])*(a-x[i]) + (b-y[i])*(b-y[i]) ));
    return tmp;
}

int main(){
    /*create a random number generator withing the range [-1,1]*/
    std::random_device rd;
    std::mt19937_64 mt(rd());
    std::uniform_real_distribution<double> dist(-1.0,1.0);
    for(unsigned int i(0);i<N;i++){
        x.push_back(dist(mt));//generate random x
        y.push_back(dist(mt));//generate random y
        /*evaluate x,y for a given function*/
        z.push_back(x.back()*y.back()+cos(x.back())*exp(y.back()*y.back()));
    }
    /*computes the matrix such that M.weights=z*/
    Matrix M(N);
    for(unsigned int i=0;i<N;i++)
        for(unsigned int j(i+1);j<N;j++)
            M(i,j) = f(sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) ));
    /*compute the invert of M. warning ! it might give wrong result if one
    line of m is 0 (all neighbours are too far) but that not where my
    problem is*/
    M.inv();
    /*using the invert of M, compute the weights*/
    for(unsigned int i(0);i<N;i++){
        weights.push_back(0.0);
        for(unsigned int j(0);j<N;j++)
            weights[i] += M(i,j)*z[j];
    }
    /*if it works, this should always return 0*/
    for(unsigned int i(0);i<N;i++)
        std::cout<<z[i]-extrapolate(x[i],y[i])<<std::endl;
}

The code is in C++ and to compile you will need the -std=c++11 -llapack flags

Answer 1

I feel so stupid... I know what is wrong...

for(unsigned int i=0;i<N;i++)
    for(unsigned int j(i);j<N;j++)
        M(i,j) = f(sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j])));

where I just changed the i+1 into i.