[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Computing distance matrix robustly
From: |
Søren Hauberg |
Subject: |
Computing distance matrix robustly |
Date: |
Sun, 02 Oct 2005 22:27:10 +0200 |
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
-------------------------------------------------------------
- Computing distance matrix robustly,
Søren Hauberg <=