[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Computing distance matrix robustly
From: |
Kevin Gross |
Subject: |
Re: Computing distance matrix robustly |
Date: |
Mon, 3 Oct 2005 09:17:16 -0400 |
søn, 02 10 2005 kl. 18:50 -0400, skrev Tom Holroyd:
On Sun, 2 Oct 2005, [utf-8] S?ren Hauberg wrote:
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.
I'm actually surprised that it worked. But perhaps (in addition
to being faster) it is true that the FPU's registers are actually
wider than a double? In that case, since X.^2 has to make
a temp array, it forces a loss of precision. I wouldn't expect
1e-7 of roundoff anyway, though. Unless the arrays in question
are floats?
I must say I was quite surprised as well. The vectorized code was:
sX2 = sum(X.^2, 2); # use sumsq instead
XX = X*X';
D = sqrt( repmat(sX2, 1, n) + repmat(sX2', m, 1) - 2*XX );
That is, I also take the sqrt. So the round off error, is only about
1e-14, before sqrt.
/Søren
Dr. Tom Holroyd
"A man of genius makes no mistakes. His errors are volitional and
are the portals of discovery." -- James Joyce
-------------------------------------------------------------
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
-------------------------------------------------------------
-------------------------------------------------------------
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
-------------------------------------------------------------