User login

Navigation

You are here

Constructing positive-definite Wendland Compact-Support RBF Matrix

Olumide's picture

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.

 

AttachmentSize
Plain text icon PolySphereConstraints.txt12.23 KB
Plain text icon cppCode.txt2.79 KB
Olumide's picture

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)

Olumide's picture

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
 

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 

Olumide's picture

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.

 

@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 

 

Olumide's picture

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

 

Subscribe to Comments for "Constructing positive-definite Wendland Compact-Support RBF Matrix"

Recent comments

More comments

Syndicate

Subscribe to Syndicate