You are here
Constructing positive-definite Wendland Compact-Support RBF Matrix
Wed, 2009-02-25 11:12 - Olumide
Further to my blog node/4883 (I wasn't aware of the forums when I started that post. Besides, the forums allow attachments)
-------------------------------------------------------
Here is the my test data(3d points) and theC++ file (rename to *.cpp) I wrote to build the matrix which is written to a file, and can be loaded into matlab or scilab as follows:
M = load("C:/tmp2/WendlandRBFMatrix.txt");
Thanks.
Attachment | Size |
---|---|
PolySphereConstraints.txt | 12.23 KB |
cppCode.txt | 2.79 KB |
Forums:
Free Tags:
Scilab version
Here is my scilab version. The resultingt matrix appears to be positive definite, so the problem must be with the matrix referencing in my C++ code.
###########################
M = load("C:/tmp2/PolySphereConstraints.txt");
n = length(M);
alpha = 1.0;
for i=1:n
x1 = M(i,1);
y1 = M(i,2);
z1 = M(i,3);
for j=1:n
x2 = M(j,1);
y2 = M(j,2);
z2 = M(j,3);
xdiff = x2 - x1;
ydiff = y2 - y1;
zdiff = z2 - z1;
r = sqrt( xdiff^2 + ydiff^2 + zdiff^2 );
rr = r/alpha;
if (rr < 1)
phi(i,j) = (1-rr)^4*(4*rr+1);
else
phi(i,j) = 0;
endif
end
end
eig(phi)
Wendland RBF Matrix
My matrix produced by my C++ code (without the constraint blocks) matches the result of my octave (not scilab) script i.e. both are positive definite. However, when I add the constraint blocks, the matrix is no longer positive definite. Below is my octave script:
#################
M = load("C:/tmp2/PolySphereConstraints.txt");
n = length(M);
alpha = 3.5;
for i=1:n
x1 = M(i,1);
y1 = M(i,2);
z1 = M(i,3);
for j=1:n
x2 = M(j,1);
y2 = M(j,2);
z2 = M(j,3);
xdiff = x2 - x1;
ydiff = y2 - y1;
zdiff = z2 - z1;
r = sqrt( xdiff^2 + ydiff^2 + zdiff^2 );
rr = r/alpha;
if (rr < 1)
matrix(i,j) = (1-rr)^4*(4*rr+1);
else
matrix(i,j) = 0;
endif
end
#constraints (refer to paper)
matrix(n+1,i) = matrix(i,n+1) = x1;
matrix(n+2,i) = matrix(i,n+2) = y1;
matrix(n+3,i) = matrix(i,n+3) = z1;
matrix(n+4,i) = matrix(i,n+4) = 1;
end
#zeros
for i=1:4
for j=1:4
matrix(n+i,n+j) = 0;
end
end
eig(matrix)
# chol(matrix) #alternatively
Re: Wendland radial basis function matrix
After glancing through the paper I can't understand what the paper means by the statement "In some cases (including the thin-plate spline solution), it is necessary
to add a first-degree polynomial P to account for the linear and constant portions of f and ensure positive definiteness of the solution".
If the matrix is equation (6) is what you're after, then the zero diagonal terms disqualify it from being positive definite (if I recall my linear algebra correctly). However, you should still be able to invert the matrix.
Could you also explain the reasoning behind the constraints:
\sum_i d_i = 0
and
\sum_i c_i d_i = 0
-- Biswajit
Re: Wendland radial basis function matrix
Those extra constraints are a relic of the variational/thin-plate spline heritage of the matrix. 'Tis a long story indeed.
BTW, I don't recall setting (all) diagonal terms to zero.
Re: Radial basis functions and disappointment
@Olumide:
First, equation (6) in the paper sets a few diagonal terms to zero. I did not say "all" :)
@ Everyone else:
This exchange is an example of what I find most disappointing about iMechanica. The original question was a request on iMechanica to debug Olumide's code. I thought the basis functions interesting and tried some of them out in Matlab and also suggested some ideas to Olumide.
After reading the paper mentioned by Olemide, I wanted to know the reason behind the choice of the constraint equations. The response was "'Tis a long story indeed'. In other words,"Please don't bother me with questions, I don't have the time to explain. But I expect people to spend the time to answer my questions". This is why iMechanica degenerates into Craig's list mode (hat tip Temesgen) from time to time until Zhigang revives it by giving it a fillip. Can't we expect better of our members?
-- Biswajit
Re: Radial basis functions and disappointment
Biswajit,
First of all I am truly grateful for your help. I really am. Your matlab code/approach helped me zero in on the problem. My reply did not imply "don't bother me with questions". I neither said so, nor was that the meaning I meant to convey. I'm sorry you read it as such.
The reason for those constraints is something I've studied for no less than 3 years. It seems most people don't really know where they come from. The scalar function phi is modelled as thin-plate spline interpolating a hight feild, and the polynomial terms were originally introduced to account for the nulll space of the bending energy of the thin-plate spline. In the stochastic framework, kriging, this is polynomial is referred to as the drift of the random function. In the variational formulation of the problem, the cosntraints ensure that the bending energy of the spline has finite values -- how exactly these constraints achieve this is the last bit I myself need to understand, and perhaps write about in detail in my dissertation. Its about time someone explained what those equations mean.
Sorry about the msunderstanding. Once again, I am truly grateful for your help.
Regards,
- Olumide
Edit: you may want to take a look at http://www-cse.ucsd.edu/classes/fa01/cse291/hhyu-presentation.pdf and http://bitsavers.org/pdf/mit/ai/aim/AIM-1430.pdf