## You are here

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

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;

}

- Olumide's blog
- Log in or register to post comments
- 5199 reads

## Comments

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

## Re: Radial basis functions

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

## http://www.math.usm.edu/csche

http://www.math.usm.edu/cschen/Matej.html