Need help/advise constructing positive-definite Wendland Compact-Support RBF Matrix

Submitted by Olumide on Tue, 02/24/2009 - 04:38

Hello,

I am trying to construct a matrix using Wendland's Compact-Support RBF as described in the following the paper, http://www.cs.umbc.edu/~rheingan/pubs/smi2001.pdf . In particular, I'm using the C2, d = 3 RBF: (1 - r)^4*(4r + 1). It seems quite simple enough, but none of the matrices I've constructed are positive definite, although Wendland's RBFs guarantee positive definite matrices -- the sparse solver I intend to use requires positive definite matrices.

Below are the relevant portions of my code in which I initialize my matrix. Is there anything I'm doing wrong?

Thanks,

- Olumide

//////////////////////////// code ////////////////////////////

for( unsigned i = 0; i < constraintCount; ++i )

{

   for( unsigned j = 0; j < i; ++j )

   {

       rbfMatrix[i][j] = rbfMatrix[j][i] = computeRBF( j , i );

   }

   rbfMatrix[i][i] = 1.0;

   rbfMatrix[constraintCount+0][i]= rbfMatrix[i][constraintCount+0] = constraint[i].x;

   rbfMatrix[constraintCount+1][i]= rbfMatrix[i][constraintCount+1] = constraint[i].y;

   rbfMatrix[constraintCount+2][i]= rbfMatrix[i][constraintCount+2] = constraint[i].z;

   rbfMatrix[constraintCount+3][i]= rbfMatrix[i][constraintCount+3] = 1.0;

}

double computeRBF( unsigned i, unsigned j )

{

   double radius2 = ( constraint[i].x - constraint[j].x )*( constraint[i].x - constraint[j].x ) + ( constraint[i].y - constraint[j].y )*( constraint[i].y - constraint[j].y ) + ( constraint[i].z - constraint[j].z )*( constraint[i].z - constraint[j].z );

   double radius = sqrt(radius2);

   double r = radius/alpha;

   return ( r < 1.0 ) ? pow( (1 - r) , 4 ) * (4*r + 1) : 0;

}

Biswajit Banerjee

Re: Radial basis functions

My Matlab version seems to work fine. See below. What's alpha in your computeRBF function? It seems to be uninitialized.

function RadialBasis

n = 10;

x = rand(1,n);

y = rand(1,n);

alpha = 1;

for i=1:length(x)

    x1 = x(i);

    y1 = y(i);

    for j=1:length(x)

      x2 = x(j);

      y2 = y(j);

      xdiff = x2-x1;

      ydiff = y2-y1;

      r = sqrt(xdiff^2 + ydiff^2);

      phi(i,j) = rbf(r, alpha);

    end

end

phi

% Check positive definite

eig(phi)

function [phi] = rbf(r, alpha)

   rr = r/alpha;

   if (rr < 1)

     phi = (1-rr)^4*(4*rr+1);

   else

     phi = 0;

   end