我已经为以下问题苦苦挣扎了两天。

我想做一个$d$一些数据的维插值。我首先尝试使用多谐波样条，但是当数据量增加时，矩阵的求逆是有问题的，因为它们几乎是奇异的（大条件数）。因此，我尝试使用紧支撑径向基函数（CSRBF）。

我无法理解的是：为什么使用完全相同的算法（见下文），它可以与薄板样条函数（或任何其他多谐波基函数）一起使用，但它与任何 CSRBF 完全失败。

在这里，我提供了一个最小的部分工作示例，它实现了两个工作功能（参见wiki）

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

和两个非工作功能（见pdf第39页）

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

这是代码，我试图尽可能地简化它，所以它效率不高，但您可以使用它并查看我的评论：

#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;
}

代码在 C++ 中，要编译，您将需要-std=c++11 -llapack标志