a-b+b != a

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

a-b+b != a

stn021
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Doug Stewart-4


On Mon, Sep 4, 2017 at 10:53 AM, stn021 <[hidden email]> wrote:
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


try this

r3 = r1+ (- r2 + r2) ;

use brackets.

--
DASCertificate for 206392


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Gordon Haverland
In reply to this post by stn021
On Mon, 4 Sep 2017 16:53:29 +0200
stn021 <[hidden email]> wrote:

> Is there a way to avoid this phenomenon ?

Learn about errors in floating point arithmetic?

This might be a start:
 https://en.wikipedia.org/wiki/Kahan_summation_algorithm

But in general you should not expect (a-b)+b to equal a.

Gord


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

James Sherman
In reply to this post by stn021
On Mon, Sep 4, 2017 at 10:53 AM, stn021 <[hidden email]> wrote:

> Hi,
>
> this is weird:
>
> r1 = round( 1e3*randn(5) ) / 1e3 ;
> r2 = round( 1e3*randn(5) ) / 1e3 ;
> r3 = r1 - r2 + r2 ;
> not_zero = r1 - r3
>
> not_zero =
>
>    0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
>    5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>    0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
>   -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
>    2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>
> So a-b+b != a
>
> The difference to zero is small, around 1e-17. But iterations can
> cause this error to increase.
>
> I use leasqr() and in each iteration the last line is
>   retval = retval - someval + someval
>
> With that additional line I get quite different results compared to
> the same program without this line even though they should be
> identical.
>
> Is there a way to avoid this phenomenon ?
>
> THX
> Stefan
>
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://lists.gnu.org/mailman/listinfo/help-octave

This is fundamental issue caused by the use of floating point numbers.
I'd suggest reading the article to get a better background in floating
point, but the key section is
https://en.wikipedia.org/wiki/Floating-point_arithmetic#Accuracy_problems

One way to avoid this issue to never test equality (or not equality)
when dealing with floating point numbers.  So instead of
>if ( a == b )
>  do stuff
>end
instead use something like:
>if ( abs(a-b) < eps(a) )
> do stuff
>end

Hope this helps,
James Sherman

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

stn021
Well, then it is easy to solve. Simply round by some power of 2:

mul = 2^20 ;
r1 = round( mul*randn(5) ) / mul ;
r2 = round( mul*randn(5) ) / mul ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0



2017-09-04 17:58 GMT+02:00 James Sherman Jr. <[hidden email]>:

> On Mon, Sep 4, 2017 at 10:53 AM, stn021 <[hidden email]> wrote:
>> Hi,
>>
>> this is weird:
>>
>> r1 = round( 1e3*randn(5) ) / 1e3 ;
>> r2 = round( 1e3*randn(5) ) / 1e3 ;
>> r3 = r1 - r2 + r2 ;
>> not_zero = r1 - r3
>>
>> not_zero =
>>
>>    0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
>>    5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>>    0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
>>   -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
>>    2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>>
>> So a-b+b != a
>>
>> The difference to zero is small, around 1e-17. But iterations can
>> cause this error to increase.
>>
>> I use leasqr() and in each iteration the last line is
>>   retval = retval - someval + someval
>>
>> With that additional line I get quite different results compared to
>> the same program without this line even though they should be
>> identical.
>>
>> Is there a way to avoid this phenomenon ?
>>
>> THX
>> Stefan
>>
>> _______________________________________________
>> Help-octave mailing list
>> [hidden email]
>> https://lists.gnu.org/mailman/listinfo/help-octave
>
> This is fundamental issue caused by the use of floating point numbers.
> I'd suggest reading the article to get a better background in floating
> point, but the key section is
> https://en.wikipedia.org/wiki/Floating-point_arithmetic#Accuracy_problems
>
> One way to avoid this issue to never test equality (or not equality)
> when dealing with floating point numbers.  So instead of
>>if ( a == b )
>>  do stuff
>>end
> instead use something like:
>>if ( abs(a-b) < eps(a) )
>> do stuff
>>end
>
> Hope this helps,
> James Sherman

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

stn021
In reply to this post by Doug Stewart-4
Hi Doug,

the idea here is to add weights to the result.

res = y + weights .* ( res-y )

If weights==1 then res should remain unchanged. Only it doesn't.
So (-r2+r2) does not really work here.

I did not include the weights in my original post because it did not seem relevant to the issue, sorry about that.

Thx, Stefan


2017-09-04 17:46 GMT+02:00 Doug Stewart <[hidden email]>:
> ...
> try this
>
> r3 = r1+ (- r2 + r2) ;
>
> use brackets.


2017-09-04 17:46 GMT+02:00 Doug Stewart <[hidden email]>:


On Mon, Sep 4, 2017 at 10:53 AM, stn021 <[hidden email]> wrote:
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


try this

r3 = r1+ (- r2 + r2) ;

use brackets.

--
DASCertificate for 206392



_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Tim Pierce
In reply to this post by stn021
No way! It's correct

in pure mathematics b - a + a = b

in computing it is not true that b -a + a = b

floating points and doubles and rounding mean that it is definitely NOT the case that mathematical identities hold

On Mon, Sep 4, 2017 at 3:53 PM, stn021 <[hidden email]> wrote:
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Mike Miller-4
In reply to this post by stn021
On Mon, Sep 04, 2017 at 18:50:29 +0200, stn021 wrote:
> the idea here is to add weights to the result.
>
> res = y + weights .* ( res-y )
>
> If weights==1 then res should remain unchanged. Only it doesn't.

If this property is important, then you should probably just add a
special case

    if (weights != 1)
      res = y + weights .* (res - y);
    endif

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

stn021
2017-09-04 19:00 GMT+02:00 Mike Miller <[hidden email]>:

> On Mon, Sep 04, 2017 at 18:50:29 +0200, stn021 wrote:
>> the idea here is to add weights to the result.
>>
>> res = y + weights .* ( res-y )
>>
>> If weights==1 then res should remain unchanged. Only it doesn't.
>
> If this property is important, then you should probably just add a
> special case
>
>     if (weights != 1)
>       res = y + weights .* (res - y);
>     endif
>
> --
> mike
Hi Mike,

not sure about this. The important property is that the weights are
applied correctly. So weights==1/2 should lead to half the value just
as weights==1 should lead to an unchanged value.

Just to be sure: weights is a vector, not a scalar. So some elements
can be 1 and others can be something else at the same time. These are
supposed to be statistical weights assigning more or less importance
to observed data.

Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Mike Miller-4
On Mon, Sep 04, 2017 at 19:18:10 +0200, stn021 wrote:
> not sure about this. The important property is that the weights are
> applied correctly. So weights==1/2 should lead to half the value just
> as weights==1 should lead to an unchanged value.
>
> Just to be sure: weights is a vector, not a scalar. So some elements
> can be 1 and others can be something else at the same time. These are
> supposed to be statistical weights assigning more or less importance
> to observed data.

Sure.

My point was that if you want to be absolutely sure that the value
doesn't change, the only solution is to avoid doing any arithmetic on
it.

It's up to you to decide how much you care about the value being exactly
the same or close enough to the same when this operation is done with a
weight of 1.

If I had to guess, res + 1e-16 is close enough to res to be considered
equal for your application.

It's not clear to me if you are asking why this limitation in precision
occurs in your expression, or if you know about floating point
arithmetic and are asking for a workaround to ensure that the value
doesn't change, or something else.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Tim Pierce
In reply to this post by Tim Pierce
TL DR; addition, subtraction, mult and div do not hold with mathematical identities

Addition, subtraction, mult and div all run the risk of 'losing accuracy' 

So every one of those operations can 'break' an identity contrasted with pure maths

Put another way, 

Try this to illustrate to yourself

        > 1e-20 + 1e20

or 

        > 2e-20 + 2e20 - 2e20

or if you really want to 'blow your mind'

        > -0

        > 0

This is worth knowing, so if you need to perform  

        ( x /10000 )* 8

 for example it will 9 times out of 10 be more 'accurate' to compute

        ( x * 8 ) / 10000

because div by a large amount is more likely to 'truncate' detail than mult ...




On Mon, Sep 4, 2017 at 5:56 PM, Tim Pierce <[hidden email]> wrote:
No way! It's correct

in pure mathematics b - a + a = b

in computing it is not true that b -a + a = b

floating points and doubles and rounding mean that it is definitely NOT the case that mathematical identities hold

On Mon, Sep 4, 2017 at 3:53 PM, stn021 <[hidden email]> wrote:
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave



_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

Tim Pierce
Sorry, I'm taking this too far and off topic, but 'interesting result' :

    >    -0 ^ -0

Anyone guess the (octave) answer without running it ? :)

On Mon, Sep 4, 2017 at 8:41 PM, Tim Pierce <[hidden email]> wrote:
TL DR; addition, subtraction, mult and div do not hold with mathematical identities

Addition, subtraction, mult and div all run the risk of 'losing accuracy' 

So every one of those operations can 'break' an identity contrasted with pure maths

Put another way, 

Try this to illustrate to yourself

        > 1e-20 + 1e20

or 

        > 2e-20 + 2e20 - 2e20

or if you really want to 'blow your mind'

        > -0

        > 0

This is worth knowing, so if you need to perform  

        ( x /10000 )* 8

 for example it will 9 times out of 10 be more 'accurate' to compute

        ( x * 8 ) / 10000

because div by a large amount is more likely to 'truncate' detail than mult ...




On Mon, Sep 4, 2017 at 5:56 PM, Tim Pierce <[hidden email]> wrote:
No way! It's correct

in pure mathematics b - a + a = b

in computing it is not true that b -a + a = b

floating points and doubles and rounding mean that it is definitely NOT the case that mathematical identities hold

On Mon, Sep 4, 2017 at 3:53 PM, stn021 <[hidden email]> wrote:
Hi,

this is weird:

r1 = round( 1e3*randn(5) ) / 1e3 ;
r2 = round( 1e3*randn(5) ) / 1e3 ;
r3 = r1 - r2 + r2 ;
not_zero = r1 - r3

not_zero =

   0.0000e+00   0.0000e+00   5.5511e-17   0.0000e+00   1.1102e-16
   5.5511e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -2.0817e-16
  -5.5511e-17   0.0000e+00   0.0000e+00  -1.1102e-16   0.0000e+00
   2.7756e-17   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00

So a-b+b != a

The difference to zero is small, around 1e-17. But iterations can
cause this error to increase.

I use leasqr() and in each iteration the last line is
  retval = retval - someval + someval

With that additional line I get quite different results compared to
the same program without this line even though they should be
identical.

Is there a way to avoid this phenomenon ?

THX
Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave




_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: a-b+b != a

stn021
Hi all,


thank you for all the input. I guess it might be a good idea to
explain in more detail what the question was about in the first place.


I want to match observed dependent and independent variables with a
non-linear model. So far I found it helpful to test a new version of
the model against a recent one in order to see if it performs better.


For that I create a new version of the program and add some new
feature to it. Then the first step is to set the parameters of the new
model so that it does the same as the previous one and the check if
the results are exactly the same. That way I can check that I did not
mess anything up.

In this case I had the idea to add variable weights that can change
from one iteration to the next. (Normally weights are fixed at the
start of the program and do not change later.)


The obvious first check of the correctness of the code is to set all
the weights to value "one". Then the results should be exactly
identical to the previous version. Only then does it make sense to
change the weights and see if the new model perfoms better than the
previous one.

Comparison gets difficult if the new program produces seemingly
randomly different results because small errors accumulate
differently. Then I cannot be sure if a slightly "better" result is
just coincidence of if the new model is really better.


The errors for a single iteration were in the range around 1e-16. So I
quickly got the idea that the issue of the different results could be
resolved by rounding, say 9 or 12 digits. But that only led to yet
another set of results.

I had forgotten to consider the fact that floating point arithmetic is
basically binary, not decimal. So rounding with 1e+9 or 1e+12 did not
really help. But rounding with 2^20 or more works just fine. Both
models give the same result.

Now I round with 2^45 ( newval = round(2^45*oldval)/2^45 ), that is
approximately 3.5e13 or 13.5 decimal digits. Any more and the
different results are back.

Stefan

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave