Re: Computing distance matrix robustly

[Kevin Gross]

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