FYI: sum improved

[Jaroslav Hajek-2]

the built-in "sum" function was improved:
http://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1

"sum" now accepts an 'extra' option that makes it use a special
summation algorithm, compensating for loss of precision (slightly more
sophisticated than the well-known Kahan's). The effect is demonstrable
by the following script:

ans =  0.015097
Elapsed time is 0.000108957 seconds.
ans = 0
Elapsed time is 0.0004 seconds.
ans =  10.015
Elapsed time is 2.69e-05 seconds.
ans =  10
Elapsed time is 3.41e-05 seconds.
ans = -7.1538
Elapsed time is 0.000391 seconds.
ans =  1.1546e-14
Elapsed time is 0.0005829 seconds.
ans =  2.8462
Elapsed time is 0.0002379 seconds.
ans =  10.000
Elapsed time is 0.000573 seconds.
ans =  24.671
Elapsed time is 0.0269 seconds.
ans = -1.6911e-12
Elapsed time is 0.06166 seconds.
ans =  34.671
Elapsed time is 0.0352 seconds.
ans = 10.0000
Elapsed time is 0.06382 seconds.

as you can see, the compensated summation always gives a much more
precise result (should be 0 or 10), at the cost of speed penalty
factor up to some 2+ (this depends on memory management). This is with
SSE arithmetics on 64-bit machines; I'm wondering what 387 arithmetics
would do. The "extra" option is not yet working for sparse matrices.

the double-accumulated non-native reduction (default for integers,
specified by "double") is now performed inline, without temporary
conversion to double if possible. This brings it on par with native
reductions, as demonstrated by the following benchmark (set n to
suitable value):

n = 2500;
disp ("single");
a = single (rand (n));
tic; sum (a, 1, "native"); toc;
tic; sum (a, 1, "double"); toc;
tic; sum (a, 2, "native"); toc
tic; sum (a, 2, "double"); toc
tic; sum (a(:), "native"); toc
tic; sum (a(:), "double"); toc
disp ("single complex");
a = complex (single (rand (n)), single (rand (n)));
tic; sum (a, 1, "native"); toc;
tic; sum (a, 1, "double"); toc;
tic; sum (a, 2, "native"); toc
tic; sum (a, 2, "double"); toc
tic; sum (a(:), "native"); toc
tic; sum (a(:), "double"); toc
disp ("int32");
a = int32 (1e4 * (rand (n) - 0.5));
tic; sum (a, 1, "native"); toc;
tic; sum (a, 1, "double"); toc;
tic; sum (a, 2, "native"); toc
tic; sum (a, 2, "double"); toc
tic; sum (a(:), "native"); toc
tic; sum (a(:), "double"); toc
disp ("uint16");
a = uint16 (1e1 * rand (n));
tic; sum (a, 1, "native"); toc;
tic; sum (a, 1, "double"); toc;
tic; sum (a, 2, "native"); toc
tic; sum (a, 2, "double"); toc
tic; sum (a(:), "native"); toc
tic; sum (a(:), "double"); toc
disp ("int64");
a = int64 (1e8 * (rand (n) - 0.5));
tic; sum (a, 1, "native"); toc;
tic; sum (a, 1, "double"); toc;
tic; sum (a, 2, "native"); toc
tic; sum (a, 2, "double"); toc
tic; sum (a(:), "native"); toc
tic; sum (a(:), "double"); toc

with a recent tip, I get (Core 2 Duo @ 2.83 GHz, g++ -O3 -march=native)

single
Elapsed time is 0.00736904 seconds.
Elapsed time is 0.0398569 seconds.
Elapsed time is 0.00691319 seconds.
Elapsed time is 0.0457029 seconds.
Elapsed time is 0.00729012 seconds.
Elapsed time is 0.0429699 seconds.
single complex
Elapsed time is 0.009444 seconds.
Elapsed time is 0.122845 seconds.
Elapsed time is 0.0125308 seconds.
Elapsed time is 0.124614 seconds.
Elapsed time is 0.00960779 seconds.
Elapsed time is 0.125515 seconds.
int32
Elapsed time is 0.00868487 seconds.
Elapsed time is 0.04076 seconds.
Elapsed time is 0.00867987 seconds.
Elapsed time is 0.0458021 seconds.
Elapsed time is 0.00956798 seconds.
Elapsed time is 0.0463171 seconds.
uint16
Elapsed time is 0.00492692 seconds.
Elapsed time is 0.039608 seconds.
Elapsed time is 0.00606394 seconds.
Elapsed time is 0.0457959 seconds.
Elapsed time is 0.00865698 seconds.
Elapsed time is 0.0426869 seconds.
int64
Elapsed time is 0.010695 seconds.
Elapsed time is 0.0449679 seconds.
Elapsed time is 0.0123532 seconds.
Elapsed time is 0.0477431 seconds.
Elapsed time is 0.0102439 seconds.
Elapsed time is 0.04808 seconds.

with the new patch, I get:

single
Elapsed time is 0.00735211 seconds.
Elapsed time is 0.00744104 seconds.
Elapsed time is 0.0052731 seconds.
Elapsed time is 0.0054481 seconds.
Elapsed time is 0.0071578 seconds.
Elapsed time is 0.00739503 seconds.
single complex
Elapsed time is 0.0096209 seconds.
Elapsed time is 0.0110161 seconds.
Elapsed time is 0.0124722 seconds.
Elapsed time is 0.0198781 seconds.
Elapsed time is 0.010736 seconds.
Elapsed time is 0.0117328 seconds.
int32
Elapsed time is 0.00881386 seconds.
Elapsed time is 0.00739193 seconds.
Elapsed time is 0.00970101 seconds.
Elapsed time is 0.0179269 seconds.
Elapsed time is 0.00954795 seconds.
Elapsed time is 0.00770593 seconds.
uint16
Elapsed time is 0.00493503 seconds.
Elapsed time is 0.00696611 seconds.
Elapsed time is 0.00503612 seconds.
Elapsed time is 0.013629 seconds.
Elapsed time is 0.00915194 seconds.
Elapsed time is 0.00722909 seconds.
int64
Elapsed time is 0.010839 seconds.
Elapsed time is 0.00875592 seconds.
Elapsed time is 0.0108941 seconds.
Elapsed time is 0.015064 seconds.
Elapsed time is 0.0130761 seconds.
Elapsed time is 0.010613 seconds.

as you can see, some operations' speed improved by a factor of 5 or more.

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[Søren Hauberg]

tir, 13 10 2009 kl. 12:38 +0200, skrev Jaroslav Hajek:
> the built-in "sum" function was improved:
> http://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1

Nice :-) It's always great to see improvements to the basic stuff.

> "sum" now accepts an 'extra' option that makes it use a special
> summation algorithm, compensating for loss of precision (slightly more
> sophisticated than the well-known Kahan's).

The name of the 'extra' option seems somewhat odd to me. Perhaps it
should be called 'extra-precise', 'improve-accuracy', or something like
that?

Søren

[Jaroslav Hajek-2]

On Tue, Oct 13, 2009 at 1:13 PM, Søren Hauberg <soren@...> wrote:

> tir, 13 10 2009 kl. 12:38 +0200, skrev Jaroslav Hajek:
>> the built-in "sum" function was improved:
>> http://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1
>
> Nice :-) It's always great to see improvements to the basic stuff.
>
>> "sum" now accepts an 'extra' option that makes it use a special
>> summation algorithm, compensating for loss of precision (slightly more
>> sophisticated than the well-known Kahan's).
>
> The name of the 'extra' option seems somewhat odd to me. Perhaps it
> should be called 'extra-precise', 'improve-accuracy', or something like
> that?
>

Well, since the option argument is currently only used to indicate
summation precision ("native" or "double"), I thought "extra" was
pretty good. I'm not against changing, but I would prefer a single
word if possible ("compensated"?). Multiword options should probably
use CamelCasing, like they do in many other functions.

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[John W. Eaton-3]

On 14-Oct-2009, Jaroslav Hajek wrote:

| Multiword options should probably
| use CamelCasing, like they do in many other functions.

Please don't force StudlyCaps on Octave users.  That is only done for
compatibility with Matlab, and in all cases that I know of, Matlab
also accepts notcamelcasing (and probably SCREAMINGNOTCAMELCASING as
well).  I'd also be happy to accept camel-casing as an option
(simply drop the hyphens, then do the comparison), but it would be
more work, and might cause some confusion among users, so I'm not sure
we should do it.

jwe

[Søren Hauberg]

ons, 14 10 2009 kl. 07:52 +0200, skrev Jaroslav Hajek:

> On Tue, Oct 13, 2009 at 1:13 PM, Søren Hauberg <soren@...> wrote:
> > tir, 13 10 2009 kl. 12:38 +0200, skrev Jaroslav Hajek:
> >> the built-in "sum" function was improved:
> >> http://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1
> >
> > Nice :-) It's always great to see improvements to the basic stuff.
> >
> >> "sum" now accepts an 'extra' option that makes it use a special
> >> summation algorithm, compensating for loss of precision (slightly more
> >> sophisticated than the well-known Kahan's).
> >
> > The name of the 'extra' option seems somewhat odd to me. Perhaps it
> > should be called 'extra-precise', 'improve-accuracy', or something like
> > that?
> >
>
> Well, since the option argument is currently only used to indicate
> summation precision ("native" or "double"), I thought "extra" was
> pretty good. I'm not against changing, but I would prefer a single
> word if possible ("compensated"?). Multiword options should probably
> use CamelCasing, like they do in many other functions.

Okay, I see the reason behind calling it 'extra'. My objection is simply
that I interpreted it is 'extra-double' or 'extra-native', which didn't
make that much sense to me. Since the word 'accurate' is used both in
the help text and in the title of the paper you followed, perhaps
'accurate' would be a good name of the option. 'compensated' also sounds
fine to me.

Søren

[Michael D. Godfrey]

On 10/14/2009 11:37 PM, Søren Hauberg wrote:
> "sum" now accepts an 'extra' option that makes it use a special
>>  >>  summation algorithm, compensating for loss of precision (slightly more
>>  >>  sophisticated than the well-known Kahan's).

Since Jaroslav did the work, I think he should decide
on the name for the extra precision.

So, lets get it done.

Michael

[Jaroslav Hajek-2]

On Thu, Oct 15, 2009 at 8:37 AM, Søren Hauberg <soren@...> wrote:

> ons, 14 10 2009 kl. 07:52 +0200, skrev Jaroslav Hajek:
>> On Tue, Oct 13, 2009 at 1:13 PM, Søren Hauberg <soren@...> wrote:
>> > tir, 13 10 2009 kl. 12:38 +0200, skrev Jaroslav Hajek:
>> >> the built-in "sum" function was improved:
>> >> http://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1
>> >
>> > Nice :-) It's always great to see improvements to the basic stuff.
>> >
>> >> "sum" now accepts an 'extra' option that makes it use a special
>> >> summation algorithm, compensating for loss of precision (slightly more
>> >> sophisticated than the well-known Kahan's).
>> >
>> > The name of the 'extra' option seems somewhat odd to me. Perhaps it
>> > should be called 'extra-precise', 'improve-accuracy', or something like
>> > that?
>> >
>>
>> Well, since the option argument is currently only used to indicate
>> summation precision ("native" or "double"), I thought "extra" was
>> pretty good. I'm not against changing, but I would prefer a single
>> word if possible ("compensated"?). Multiword options should probably
>> use CamelCasing, like they do in many other functions.
>
> Okay, I see the reason behind calling it 'extra'. My objection is simply
> that I interpreted it is 'extra-double' or 'extra-native', which didn't
> make that much sense to me.

No, "extra" can't be specified together with "double" or "native",
it's just a third option for the working precision. Maybe it could be
clarified in the docstring. My idea was that you simply specify that
you want the sum to be extra precise. For single precision, this means
it's computed in double (same as 'double'), for double precision, the
compensated algorithm is used.

> Since the word 'accurate' is used both in
> the help text and in the title of the paper you followed, perhaps
> 'accurate' would be a good name of the option. 'compensated' also sounds
> fine to me.
>

To me, "accurate" doesn't seem any better than "extra". "compensated"
is pretty unambiguous, but slightly misleading, because using "extra"
on single prec numbers doesn't use the compensated summation in single
precision (that would make little sense, because summing in double is
better). I'm still unconvinced to abandon "extra", but I can easily
make alternative spellings like "extraPrecision" or "accurate".
I'll try to clarify the docstring a bit, though.

regards

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz