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 |
On Mon, Sep 4, 2017 at 10:53 AM, stn021 <[hidden email]> wrote: Hi, try this r3 = r1+ (- r2 + r2) ; use brackets. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
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 |
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 |
In reply to this post by Doug Stewart-4
Hi Doug, I did not include the weights in my original post because it did not seem relevant to the issue, sorry about that.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. 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]>:
_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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, _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
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 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 |
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 |
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:
_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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:
_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
Free forum by Nabble | Edit this page |