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