Re: Computing distance matrix robustly

[Søren Hauberg]

Hi again,
Sorry about replying to my own post, but Tom Holroyd, contacted me
off-list with a solution. Simply replace

sX2 = sum(X.^2, 2);

with

sX2 = sumsq(X, 2);

in the vectorized version, and then things work.

/Søren

søn, 02 10 2005 kl. 22:27 +0200, skrev Søren Hauberg:
> Hi,
> This is a long post, so feel free to spik it.
> I'm trying to compute this distance matrix of a set of points. That is,
> I have an MxD matrix consisting of M points in D dimensional space, and
> I want to compute the distances between all the points. This distance
> matrix can then be defined as,
> 
> D(i,j) = norm( X(i,:) - X(j,:) )
> 
> where X is the above mentioned matrix. This can be implemented using two
> for-loops, but this is very slow, so I want to vectorize it. We note
> that,
> 
> D(i,j) = sqrt(sum( (X(i,:) - X(j,:)).^2 ))
>        = sqrt(sum( X(i,:).^2 + X(j,:).^2 - 2*X(i,j)^2 ))
> 
> A vectorized version of the computation then becomes,
> 
> sX2 = sum(X.^2, 2);
> XX = X*X';
> D = sqrt( repmat(sX2, 1, n) + repmat(sX2', m, 1) - 2*XX );
> 
> Now this is fast, but it's not robust. The diagonal element of this
> matrix should be zero (the distance between a point and it self is
> zero), but they sometimes get values that are approx. 4*10^-7 and
> sometimes they are imaginary. Now I could force the diagonal elements to
> be zero, but I would like the function to be more general, so that it
> works with two matrices instead of one, which might result in a
> non-square distance matrix, so in general that's not an option.
> 
> Does anybody have an idea on how to fix this?
> 
> Thanks,
> Søren
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------
> 



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------