Thread Subject: Dumbfounding Loop Result

Hi everyone -

I am working on a loop and came across a very confusing
result, which is captured in the code below (my code is
more complex, but this captures the problem cleanly).

When I add q+r+s and compare it to s+q+r, I get a different
result for about 20% of values of q,r,s.

Is this a numerical approximation issue? If so, is there
an "almost equal to" equality operator or something?

Any help would be greatly appreciated...I am confused.

>>>>>

% code returns 'incorrect' many times

for q=0:.1:1
    for r=0:.1:1
        for s=0:.1:1
            
            sum1=q+r+s;
            sum2=s+q+r;
            
            if sum1==sum2
                correct=1;
            else
                incorrect=1
            end

        end
    end
end

>>>>>

"Ned A" <ned456@hotmail.com> wrote in message
<g4e1pb$2qt$1@fred.mathworks.com>...
> Hi everyone -
>
> I am working on a loop and came across a very confusing
> result, which is captured in the code below (my code is
> more complex, but this captures the problem cleanly).
>
> When I add q+r+s and compare it to s+q+r, I get a
different
> result for about 20% of values of q,r,s.
>
> Is this a numerical approximation issue? If so, is there
> an "almost equal to" equality operator or something?
>
> Any help would be greatly appreciated...I am confused.
>
> >>>>>
>
> % code returns 'incorrect' many times
>
> for q=0:.1:1
> for r=0:.1:1
> for s=0:.1:1
>
> sum1=q+r+s;
> sum2=s+q+r;
>
> if sum1==sum2
> correct=1;
> else
> incorrect=1
> end
>
> end
> end
> end
>
> >>>>>
>

It's probably a numerical accuracy issue. You should try
something like 'abs(sum1-sum2) < eps' to check for equality.

Ken

"Ned A" <ned456@hotmail.com> wrote in message
<g4e1pb$2qt$1@fred.mathworks.com>...
< ... SNIP fp issues

This is typical behavior of floating point arithmetic.
Search this newsgroup for numerous other examples, such as

.1 + .2 + .3 - .6

It arises from the fact that finite computer systems cannot
store numbers like .1 exactly, and therefore introduce
small round-off errors:

fprintf('%.20f',.6)

to check for equality, use something like:

abs(fp1 - fp2) < threshold

hth
Jos

In article <g4e1pb$2qt$1@fred.mathworks.com>, Ned A <ned456@hotmail.com> wrote:

>I am working on a loop and came across a very confusing
>result, which is captured in the code below (my code is
>more complex, but this captures the problem cleanly).

>When I add q+r+s and compare it to s+q+r, I get a different
>result for about 20% of values of q,r,s.

>Is this a numerical approximation issue?

Not exactly: the issue is that '+' is not a commutative operation
with any floating point system that uses a fixed-length mantissa
togethr with a variable exponent. The same problem occurs
with processors with decimal floating point arithmetic
(though perhaps not with those specific loop values.)


>If so, is there
>an "almost equal to" equality operator or something?

No there isn't.

I am not an expert on floating point arithmetic, so the
following might not be exactly correct, but the "approximately equal"
operator is close to:

  abs(a - b) < eps(max(abs(a),abs(b)))

If you know that you do not require high precision, then you can
approximate this with (e.g.)

  abs(a-b) < 1e-6

or whatever tolerance level is appropriate for the problem.


Note:
  abs(a-b) < eps
is -not- correct, not unless your tolerance just happens to be
the fixed number eps(1) ( about 2.22045e-16 ). Using eps
instead of eps() of the expression fails to take into account
*relative* tolerances, such as the fact that beyond 2^51 a
single bit difference in a double precision number results in
a change of more than 1 in the result.
--
  "The beauties of conception are always superior to those of
   expression." -- Walter J. Phillips

> I am not an expert on floating point arithmetic, so the
> following might not be exactly correct, but
the "approximately equal"
> operator is close to:
>
> abs(a - b) < eps(max(abs(a),abs(b)))
>
> If you know that you do not require high precision, then
you can
> approximate this with (e.g.)
>
> abs(a-b) < 1e-6
>
> or whatever tolerance level is appropriate for the
problem.
>
>
> Note:
> abs(a-b) < eps
> is -not- correct, not unless your tolerance just happens
to be
> the fixed number eps(1) ( about 2.22045e-16 ). Using eps
> instead of eps() of the expression fails to take into
account
> *relative* tolerances, such as the fact that beyond 2^51 a
> single bit difference in a double precision number
results in
> a change of more than 1 in the result.

Good point, Walter. I forgot to add that extra piece to the
equation.

Ken

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in message <g4e2vc
$cnt$1@canopus.cc.umanitoba.ca>...
> ......
> Not exactly: the issue is that '+' is not a commutative operation
> with any floating point system that uses a fixed-length mantissa
> togethr with a variable exponent. The same problem occurs
> with processors with decimal floating point arithmetic
> (though perhaps not with those specific loop values.)
> ......

  If the floating point system used adheres to the IEEE 754 standard, addition
and multiplication always satisfy the commutative law. Binary operations a+b
and b+a will always give exactly equal results, and the same is true of a*b
and b*a. This is carefully spelled out in the standard. It is the associative law
involving three successive operands that is not strictly obeyed. The result of
(a+b)+c is not always exactly equal to a+(b+c). This is because different
intermediate results are being rounded in the two different groupings.

  Also it should be pointed out that the above is not the result of inexact
representations of decimal fractions such as .1, .2, etc., (which are
responsible for an entirely different class of errors.) Ned's two operations q
+r+s and s+q+r will give exactly the same values as (r+q)+s and r+(q+s),
respectively, (courtesy of the commutative law.) However these latter two can
differ because the associative law is not always followed precisely, and this
has nothing to do with whether q, r, and s may be inexact approximations of
decimal fractions.

Roger Stafford

In article <g4ea5d$6oc$1@fred.mathworks.com>,
Roger Stafford <ellieandrogerxyzzy@mindspring.com.invalid> wrote:
>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in message <g4e2vc
>$cnt$1@canopus.cc.umanitoba.ca>...

>> Not exactly: the issue is that '+' is not a commutative operation
>> with any floating point system that uses a fixed-length mantissa
>> togethr with a variable exponent.

>It is the associative law
>involving three successive operands that is not strictly obeyed.

You are correct, I got the name of the law confused.
--
  "What we have to do is to be forever curiously testing new
  opinions and courting new impressions." -- Walter Pater

"Ned A" <ned456@hotmail.com> wrote in message <g4e1pb$2qt
$1@fred.mathworks.com>...
> .......
> When I add q+r+s and compare it to s+q+r, I get a different
> result for about 20% of values of q,r,s.
> ........

  Ned, I stated that matlab could get different results from (a+b)+c than from
a+(b+c), but did not get around to giving an example to show how this can
occur. It is perhaps easiest to see with integers. Suppose you define the
three very large integers:

 a = 7000000000000002;
 b = 6000000000000001;
 c = -5000000000000000;
 
The quantities a, b, and (absolute value of) c all lie in the range from 2^52 to
2^53-1 where all numbers in matlab's double precision floating point
('double') must in fact be integers. If we first compute a+b, we will get a
larger number lying in the range 2^53 to 2^54-1 where all numbers must be
even integers. (This is because it would require 54 bits to express an odd
integer in this range, and 'double' has only 53 bits available for this purpose.)
Since the true sum of a and b, 13000000000000003, is odd, this sum must
therefore be rounded to an even value to give a valid 'double' number. It is
rounded up by 1 (to the nearest multiple of 4):

 a+b = 13000000000000004

Then adding c to this gives

 (a+b)+c = 8000000000000004

which is 1 higher than the true answer. On the other hand, if we compute the
sum of b and c first and then add a, all quantities remain in the range, 2^52
to 2^53-1 where no rounding adjustment is necessary and the exactly correct
(and different) answer

 a+(b+c) = 8000000000000003

is obtained.

  This demonstrates how the associative law of addition can fail with floating
point numbers. Moreover it is a phenomenon that, as you have learned to
your dismay, carries over to all ranges of floating point numbers, integer and
fractional. Out of any three numbers, a, b, and c, calculating a+b first and
then adding c produces the intermediate quantity a+b, while adding b and c
first and then adding a, yields the intermediate quantity b+c which may be
very different from a+b and therefore possibly possess a different rounding
situation, as demonstrated in the above example.

  As you can perhaps see from this example, this is an inherent characteristic
of all floating point systems, binary, decimal, or otherwise, and not a
peculiarity of matlab or the IEEE 754 standard. There is simply no practical
way of avoiding this kind of rounding discrepancy for addition in any floating
point system having a fixed number of bits.

  A similar situation also holds for the product of three numbers - the
associative law of multiplication is often violated if you are looking for exact
equality, but of course only to the extent of about one part in 2^52 or so
parts.

Roger Stafford

"Roger Stafford"
<ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g4hk75$jd$1@fred.mathworks.com>...
> "Ned A" <ned456@hotmail.com> wrote in message <g4e1pb$2qt
> $1@fred.mathworks.com>...
> > .......
> > When I add q+r+s and compare it to s+q+r, I get a
different
> > result for about 20% of values of q,r,s.
> > ........
>
> Ned, I stated that matlab could get different results
from (a+b)+c than from
> a+(b+c), but did not get around to giving an example to
show how this can
> occur. It is perhaps easiest to see with integers.
Suppose you define the
> three very large integers:
>
> a = 7000000000000002;
> b = 6000000000000001;
> c = -5000000000000000;
>
> The quantities a, b, and (absolute value of) c all lie in
the range from 2^52 to
> 2^53-1 where all numbers in matlab's double precision
floating point
> ('double') must in fact be integers. If we first compute
a+b, we will get a
> larger number lying in the range 2^53 to 2^54-1 where all
numbers must be
> even integers. (This is because it would require 54 bits
to express an odd
> integer in this range, and 'double' has only 53 bits
available for this purpose.)
> Since the true sum of a and b, 13000000000000003, is odd,
this sum must
> therefore be rounded to an even value to give a
valid 'double' number. It is
> rounded up by 1 (to the nearest multiple of 4):
>
> a+b = 13000000000000004
>
> Then adding c to this gives
>
> (a+b)+c = 8000000000000004
>
> which is 1 higher than the true answer. On the other
hand, if we compute the
> sum of b and c first and then add a, all quantities
remain in the range, 2^52
> to 2^53-1 where no rounding adjustment is necessary and
the exactly correct
> (and different) answer
>
> a+(b+c) = 8000000000000003
>
> is obtained.
>
> This demonstrates how the associative law of addition
can fail with floating
> point numbers. Moreover it is a phenomenon that, as you
have learned to
> your dismay, carries over to all ranges of floating point
numbers, integer and
> fractional. Out of any three numbers, a, b, and c,
calculating a+b first and
> then adding c produces the intermediate quantity a+b,
while adding b and c
> first and then adding a, yields the intermediate quantity
b+c which may be
> very different from a+b and therefore possibly possess a
different rounding
> situation, as demonstrated in the above example.
>
> As you can perhaps see from this example, this is an
inherent characteristic
> of all floating point systems, binary, decimal, or
otherwise, and not a
> peculiarity of matlab or the IEEE 754 standard. There is
simply no practical
> way of avoiding this kind of rounding discrepancy for
addition in any floating
> point system having a fixed number of bits.
>
> A similar situation also holds for the product of three
numbers - the
> associative law of multiplication is often violated if
you are looking for exact
> equality, but of course only to the extent of about one
part in 2^52 or so
> parts.


The Fortran die-hards often use this as one argument to
further their cause. Quoting from a rather
amusing "COMPARISON OF FORTRAN AND C" web page:

      1) The order of evaluation of arithmetical expressions
          is defined precisely, and can be controlled with
          parentheses.

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <g4hk75$jd$1@fred.mathworks.com>...
[snip]
> Ned, I stated that matlab could get different results
from (a+b)+c than from
> a+(b+c), but did not get around to giving an example to
show how this can
> occur. It is perhaps easiest to see with integers.
Suppose you define the
> three very large integers:
>
> a = 7000000000000002;
> b = 6000000000000001;
> c = -5000000000000000;
>
> The quantities a, b, and (absolute value of) c all lie in
the range from 2^52 to
> 2^53-1 where all numbers in matlab's double precision
floating point
> ('double') must in fact be integers. If we first compute
a+b, we will get a
> larger number lying in the range 2^53 to 2^54-1 where all
numbers must be
> even integers. (This is because it would require 54 bits
to express an odd
> integer in this range, and 'double' has only 53 bits
available for this purpose.)
> Since the true sum of a and b, 13000000000000003, is odd,
this sum must
> therefore be rounded to an even value to give a valid
'double' number. It is
> rounded up by 1 (to the nearest multiple of 4):
>
> a+b = 13000000000000004
>
> Then adding c to this gives
>
> (a+b)+c = 8000000000000004
>
> which is 1 higher than the true answer. On the other
hand, if we compute the
> sum of b and c first and then add a, all quantities remain
in the range, 2^52
> to 2^53-1 where no rounding adjustment is necessary and
the exactly correct
> (and different) answer
>
> a+(b+c) = 8000000000000003
>
> is obtained.
>
> This demonstrates how the associative law of addition
can fail with floating
> point numbers. Moreover it is a phenomenon that, as you
have learned to
> your dismay, carries over to all ranges of floating point
numbers, integer and
> fractional. Out of any three numbers, a, b, and c,
calculating a+b first and
> then adding c produces the intermediate quantity a+b,
while adding b and c
> first and then adding a, yields the intermediate quantity
b+c which may be
> very different from a+b and therefore possibly possess a
different rounding
> situation, as demonstrated in the above example.
>
> As you can perhaps see from this example, this is an
inherent characteristic
> of all floating point systems, binary, decimal, or
otherwise, and not a
> peculiarity of matlab or the IEEE 754 standard. There is
simply no practical
> way of avoiding this kind of rounding discrepancy for
addition in any floating
> point system having a fixed number of bits.
>
> A similar situation also holds for the product of three
numbers - the
> associative law of multiplication is often violated if you
are looking for exact
> equality, but of course only to the extent of about one
part in 2^52 or so
> parts.
>
> Roger Stafford

This situation gets worse when calculating variance of a
long series of data. Given a long series, x, of identical
non-zero decimal fractions (e.g. 0.1), mean(x) != 0.1, but
it is close. This estimate can be improved by
mean(x-mean(x)) + mean(x).

mean(x - mean(x)) is zero for infinite precision and it is
usually very small for double precision. There are still
pathological cases, e.g. a long series of 0.1 followed by an
equally long series of -0.1 (i.e. where the true mean is
zero), but it is usually good enough for stochastic data.

I used this technique to formulate a robust calculation of
standard deviation (and higher statistical moments). For
long series of constant values, there can still be a square
root of a very small negative number but I just replace
these with zero as it MUST have resulted from finite
precision errors.

Andy Robb.

Even with my (mean(x - mean(x)) + mean(x)) solution, there
will probably still be an error in the least significant bit
of the mantissa because of rounding in the division by n to
convert from sum() to mean() before the final addition.

The same is true if you multiply mean(x) by n to calculate a
corrected sum(x). (It is not valid to use sum(x) instead of
mean(x)*n, because we are calculating a correction to the
original estimate of mean(x).)

If you implement this in a MEX-file, just the final division
(or multiplication) and addition can be performed as 'long
double' and this will eliminate practically all differences
due to initial order - it works for all the original problem
values.