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

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

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

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

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

Olumide

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