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{
        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_];
            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);

        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*/
    /*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*/
    /*using the invert of M, compute the weights*/
    for(unsigned int i(0);i<N;i++){
        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++)

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


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.

  • 1
    $\begingroup$ There's absolutely no need to feel stupid. We all know this kind of problems. $\endgroup$ – davidhigh Jan 18 at 0:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.