User login

Navigation

You are here

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

Olumide's picture

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

Comments

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

Olumide's picture

 Thanks Biswajit , alpha is a gloal variable and it is initialized. Also, I am working in 3D. 

Please refer to my recent forum thread http://imechanica.org/node/4912 . Regards,

Olumide

Subscribe to Comments for " Need help/advise constructing positive-definite Wendland Compact-Support RBF Matrix "

Recent comments

More comments

Syndicate

Subscribe to Syndicate