Computing distance matrix robustly

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