octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: A problem with range objects and floating point numbers


From: s.jo
Subject: Re: A problem with range objects and floating point numbers
Date: Sun, 28 Sep 2014 20:20:16 -0700 (PDT)

Daniel Sebald wrote
> On 09/26/2014 02:43 AM, s.jo wrote:
>> Daniel Sebald wrote
> [snip]
>>> To be factorable doesn't necessarily require that the floating point
>>> step size be an exact number system equivalent, i.e., that .1 equal
>>> exactly 1/10.  It just requires the machine arithmetic to be as expected
>>> when producing an operation result, i.e.,
>>>
>>> octave-cli:1>  rem(-2,.1) == 0
>>> ans =  1
>>> octave-cli:2>  rem(0,.1) == 0
>>> ans =  1
>>> octave-cli:3>  [-20:1:0]*.1
>>> ans =
>>>
>>>    Columns 1 through 5:
>>>
>>>     -2.00000  -1.90000  -1.80000  -1.70000  -1.60000
>>>
>>>    Columns 6 through 10:
>>>
>>>     -1.50000  -1.40000  -1.30000  -1.20000  -1.10000
>>>
>>>    Columns 11 through 15:
>>>
>>>     -1.00000  -0.90000  -0.80000  -0.70000  -0.60000
>>>
>>>    Columns 16 through 20:
>>>
>>>     -0.50000  -0.40000  -0.30000  -0.20000  -0.10000
>>>
>>>    Column 21:
>>>
>>>      0.00000
>>>
>>>
>>> Jo, what do you get for the above scripts?
>>>
>>> Dan
>>
>> Dan and others,
>> Sorry for late reply.
>>
>> Your integerized stepping seems to be a decimation process.
>> You used step-size 1. Any bigger integers such as 10 or 10^17 are allowed
>> in
>> the decimation process.
> 
> By decimation, I take it you mean somehow converting the number to true 
> decimal representation, similar to the decimal arithmetic toolbox Oliver 
> alluded to.  I don't think that's the case.
> 
> 
>> Comparing the results of colon operators with floating-point numbers,
>> Matlab seems to use such decimation process, but Octave does not.
>> I agree that this decimation feature of colon operator is far from any
>> IEEE
>> standard.
>>
>> Also I point out that the results of linspace function and colon operator
>> in
>> Matlab are different.
>> This may imply decimation processes are different.
>> By inspection, the linspace result is much close to decimal number.
>> I guess that the decimation process of Matlab's colon operator is simpler
>> for speed.
> 
> No, I doubt that.  The linspace and colon operator are probably two 
> slightly different algorithms.  For linspace the spacing is exactly an 
> integer, by definition.  That's not necessarily so for the colon 
> operator where the spacing can result in a fraction of the step size at 
> the last interval.
> 
> Looking at the code a bit, I see that the range (colon) operator does 
> attempt to use integer multiplication as much as possible, e.g., from 
> Range::matrix_value:
> 
>        cache.resize (1, rng_nelem);
>        double b = rng_base;
>        double increment = rng_inc;
>        if (rng_nelem > 0)
>          {
>            // The first element must always be *exactly* the base.
>            // E.g, -0 would otherwise become +0 in the loop (-0 + 
> 0*increment).
>            cache(0) = b;
>            for (octave_idx_type i = 1; i < rng_nelem; i++)
>              cache(i) = b + i * increment;
>          }
> 
> So, the accumulation of addition errors doesn't seem to be an issue here.
> 
> 
>> So far, I learned that integer range can be generated 'safely' with colon
>> operator.
>> Regarding to floating-point number range, we may want to use the function
>> "linspace" for the best decimation.
>>
>> I don't know who is responsible for the colon operator source code.
>> But I think that if there is a document for comparing colon operator with
>> linspace function,
>> that will be helpful.
>>
>> There is another question: I found that (-20:1:1)*.1 and [-20:1:1]*0.1
>> have
>> different results.
>>
>>
>> The range with [...] is better decimated than the case of (...).
> 
> Interesting.  You may have it right there.  And in fact, it appears the 
> conversation in the bug report
> 
> https://savannah.gnu.org/bugs/?42589
> 
> turns to this issue.  John describes how a "Range" class retains the 
> description of the range rather than converting it to a matrix class. 
> That actually may be a good idea because it could possibly allow 
> retaining double precision if needed somehow (e.g., machines where the 
> single precision is markedly different from "double").  Whereas 
> converting the range to a matrix might make the result single precision 
> depending on compiler settings or whatnot.  I'm just speaking in the 
> abstract without any example.
> 
> So, in your example, the difference is (-20:1:1) is a Range description 
> (base, increment, limit), where [-20:1:1] is a matrix_value 
> implementation of that range (all the actual 21 values).
> 
> Let's try to discern which of these might be the more accurate result:
> 
>>> t_1 = (-20:1:1)*.1;
>>> t_2 = [-20:1:1]*.1;
>>> t_f = [-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 ...
>         -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1];
>>> t_d = double(zeros(1,21));
>>> t_d(1) = -2; t_d(2) = -1.9; t_d(3) = -1.8; t_d(4) = -1.7; t_d(5) = -1.6;
>>> t_d(6) = -1.5; t_d(7) = -1.4; t_d(8) = -1.3; t_d(9) = -1.2; t_d(10) =
>>> -1.1;
>>> t_d(11) = -1; t_d(12) = -0.9; t_d(13) = -0.8; t_d(14) = -0.7; t_d(15) =
>>> -0.6;
>>> t_d(16) = -0.5; t_d(17) = -0.4; t_d(18) = -0.3; t_d(19) = -0.2; t_d(20)
>>> = -0.1;
>>> t_d(21) = 0; t_d(22) = 0.1;
>>> max(abs(t_d - t_f))
> ans = 0
>>> norm(abs(t_d - t_f))
> ans = 0
>>> max(abs(t_d - t_1))
> ans =    2.2204e-16
>>> norm(abs(t_d - t_1))
> ans =    4.3088e-16
>>> max(abs(t_d - t_2))
> ans =    2.2204e-16
>>> norm(abs(t_d - t_2))
> ans =    4.7429e-16
>>> max(abs(t_1 - t_2))
> ans =    2.2204e-16
>>> [abs(t_d-t_1)' abs(t_d-t_2)']
> ans =
> 
>     0.0000e+00   0.0000e+00
>     0.0000e+00   2.2204e-16
>     0.0000e+00   0.0000e+00
>     0.0000e+00   2.2204e-16
>     0.0000e+00   0.0000e+00
>     0.0000e+00   0.0000e+00
>     0.0000e+00   2.2204e-16
>     2.2204e-16   0.0000e+00
>     0.0000e+00   2.2204e-16
>     0.0000e+00   0.0000e+00
>     0.0000e+00   0.0000e+00
>     1.1102e-16   0.0000e+00
>     2.2204e-16   0.0000e+00
>     0.0000e+00   1.1102e-16
>     1.1102e-16   1.1102e-16
>     0.0000e+00   0.0000e+00
>     1.1102e-16   0.0000e+00
>     1.6653e-16   5.5511e-17
>     5.5511e-17   0.0000e+00
>     1.3878e-16   0.0000e+00
>     0.0000e+00   0.0000e+00
>     0.0000e+00   0.0000e+00
> 
>  From the norm(abs(t_d - t_f)) = 0 result, I'd say my machine/compiler 
> is using double precision by default.  Also, there is a difference 
> between (A:1:B)*d and [A:1:B]*d but no implementation seems better than 
> the other.  If anything (A:1:B)*d is marginally better than [A:1:B]*d, 
> but statistically speaking the norm difference of 4.3405e-17 is probably 
> not conclusive.  However, if I were to guess, I'd say (A:1:B)*d retains 
> accuracy better because it doesn't cast to a IEEE float representation 
> and then multiply.  The computations may be all done with more bits in 
> the ALU.  But we're talking the very last bit here, I believe.
> 
> 
>> More interestingly, if I transpose the range array  such as
>> (-20:1:1)'*.1,
>> then I have the same result in [-20:1:1]'*.1.
> 
> The transpose likely forces the Range description into a matrix_value 
> implementation.
> 
> 
>> Is this a bug? Did someone say that there is a patch for this?
>> I am useing Ocatve 3.8.2 on cygwin.
> 
> OK, Cygwin.  So there's where there could be some differences in 
> compiler settings that make your computations single precision and 
> therefore showing larger errors in the arithmetic than what I'm seeing. 
>   Try the script lines above on your machine and see what you get.  You 
> probably should have the latest code with the Range patch applied.  In 
> any case, though, if you can identify a single/double precision issue, 
> it might suggest that the Octave code isn't retaining precision...or it 
> could simply be that on your system Matlab is using double precision at 
> compilation whereas Octave is not.
> 
> Dan

Dear Dan,

Thanks to your comments, I got more clear idea on ranges...
Please find my opinions as follows:

>From https://www.gnu.org/software/octave/doc/interpreter/Ranges.html , 
I can see the difference among range object, matrix object, and linspace.
I understand that the transpose function triggers converting range into
matrix type.

octave:14> (-20:1:1) - [-20:1:1]
ans =

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

octave:15> (-20:1:1)*.1 - [-20:1:1]*.1
ans =

 Columns 1 through 5:

           0  2.2204e-16           0  2.2204e-16           0

 Columns 6 through 10:

           0  2.2204e-16           0  2.2204e-16  2.2204e-16

 Columns 11 through 15:

           0  1.1102e-16  1.1102e-16  1.1102e-16  2.2204e-16

 Columns 16 through 20:

  1.1102e-16  1.1102e-16  1.1102e-16  1.1102e-16  1.1102e-16

 Columns 21 and 22:

  1.1102e-16           0


But I could not follow the above result.
It is ok that range and matrix types show the same results as they are.
Speaking of the mixed type operation with double, they look different!
There seems to be a bug in mixed type operation of range object and double
number.
I installed Octave 3.8.2 from cygwin package repository, not compiled them
from source code.
I believe that range type operated mixed with double should be same as
matrix type operated mixed with double.

Anyhow, I agree that my previous conclusion comparing (A:1:B)*d and
[A:1:B]*d was NOT correct.
Both expressions have similar performance over floating-point
representations.
There is a special script that I use when I test the ranges: 
  sum(abs(sin(k*pi))<eps), k is ranges in floating-point numbers.
This script that counts the zero-crossover points of sin is motivated from 
http://www.mathworks.co.kr/company/newsletters/articles/the-tetragamma-function-and-numerical-craftsmanship.html

>From the result of the script, I can determine that pi multiplied by decimal
points such 1.0*pi, 2.0*pi, ... 1000.0*pi are well generated from range
operations.
As a result, both (A:1:B)*d and [A:1:B]*d are poor than linspace in Octave.
However, in Matlab, all results from range type, matrix-range type and
linspace pass the 'abs(sin(k*pi))<eps' test for large k.
See the complete script and their results (comparing Octave and Matlab) as
follows:

== sin test script ==
version
% generate ranges : 
%  1. range object, 2. matrix object, 3. linspace, 4. perfect decimation
N=1000;
range_test1=(-N:0.1:N);
range_test2=[-N:0.1:N];
range_test3=linspace(-N,N,20*N+1);
range_test4=(str2num(sprintf('%.12f ',-N:0.1:N)));
% count sin(k*pi)==0 : 0 or 1 expected due to floating point representation
error
sum(sin(range_test1*pi)==0)
sum(sin(range_test2*pi)==0)
sum(sin(range_test3*pi)==0)
sum(sin(range_test4*pi)==0)
% count sin(k*pi)<eps : what we want is 2*N+1 in case of the best decimation
sum(abs(sin(range_test1*pi))<eps(range_test1*pi)) 
sum(abs(sin(range_test2*pi))<eps(range_test2*pi)) 
sum(abs(sin(range_test3*pi))<eps(range_test3*pi)) 
sum(abs(sin(range_test4*pi))<eps(range_test4*pi)) 

== Octave 3.8.2 result ==
+ version
ans = 3.8.2
+ ## generate ranges : 
+ ##  1. range object, 2. matrix object, 3. linspace, 4. perfect decimation
+ N = 1000;
+ range_test1 = (-N:0.1:N);
+ range_test2 = [-N:0.1:N];
+ range_test3 = linspace (-N, N, 20 * N + 1);
+ range_test4 = (str2num (sprintf ('%.12f ', -N:0.1:N)));
+ ## count sin(k*pi)==0 : 0 or 1 expected due to floating point
representation error
+ sum (sin (range_test1 * pi) == 0)
ans =          0
+ sum (sin (range_test2 * pi) == 0)
ans =          0
+ sum (sin (range_test3 * pi) == 0)
ans =          1
+ sum (sin (range_test4 * pi) == 0)
ans =          1
+ ## count sin(k*pi)<eps : what we want is 2*N+1 in case of the best
decimation
+ sum (abs (sin (range_test1 * pi)) < eps (range_test1 * pi))
ans =       1055
+ sum (abs (sin (range_test2 * pi)) < eps (range_test2 * pi))
ans =       1168
+ sum (abs (sin (range_test3 * pi)) < eps (range_test3 * pi))
ans =       2001
+ sum (abs (sin (range_test4 * pi)) < eps (range_test4 * pi))
ans =       2001


== Matlab 2012b result ==

version
ans =
8.0.0.783 (R2012b)
% generate ranges : 
%  1. range object, 2. matrix object, 3. linspace, 4. perfect decimation
N=1000;
range_test1=(-N:0.1:N);
range_test2=[-N:0.1:N];
range_test3=linspace(-N,N,20*N+1);
range_test4=(str2num(sprintf('%.12f ',-N:0.1:N)));
% count sin(k*pi)==0 : 0 or 1 expected due to floating point representation
error
sum(sin(range_test1*pi)==0)
ans =
     1
sum(sin(range_test2*pi)==0)
ans =
     1
sum(sin(range_test3*pi)==0)
ans =
     1
sum(sin(range_test4*pi)==0)
ans =
     1
% count sin(k*pi)<eps : what we want is 2*N+1 in case of the best decimation
sum(abs(sin(range_test1*pi))<eps(range_test1*pi)) 
ans =
        2001
sum(abs(sin(range_test2*pi))<eps(range_test2*pi)) 
ans =
        2001
sum(abs(sin(range_test3*pi))<eps(range_test3*pi)) 
ans =
        2001
sum(abs(sin(range_test4*pi))<eps(range_test4*pi)) 
ans =
        2001


== discussion ==
Range object with floating-point numbers is not much emphasized so far.

In Octave, floating-point ranges generated from range object may show poorer
performance than
those from linspace function, when evaluating trigonometric functions such
as sin, cos and exp(complex arguments).
Meanwile, Matlab show the same performance among the ranges and linspace
function 
when we apply abs(sin(k*pi))<eps test for large k. 
For example, I counted zero-crossover points in sin function in the range of
-1000*pi~1000pi,
with 0.1*pi step-size. The true value will be 2001 points in that range.
Matlab showed the correct count for any range method, but Octave showed the
correct count only in linspace function.

Note that Matlab seems to have no range type. 
Speaking of range object, there seems to be a bug with mixed type operation
with double precision numbers in Octave.

For numberical application intensively using trigonometric function with
large angles,
we had better use linspace function to generate a range of radian angles in
Octave.
----


ps. I added more test script for range hex forms and their mixed-type
operation. 
We can see how they differ from matlab's outcome.
Note that range*pi and matrix_range*pi show same in Matlab because Matlab
does not have range object.
Also note that linspace functions of Matlab and Octave are not matched even
though both pass abs(sin(k*pi))<eps test.
Octave's linspace is very close to the full decimation, but Matlab's is not. 
There seems to be an unknown additional decimal process in Matlab.

Speaking of zero-crossover check of sin function, I wrote that script to see
sin function behavior with floating-point numbers. I happen to find that the
results vary according to range methods in Octave.
Such errors occur more likely on -0 side. However, Matlab does not show such
errors.
This drove me here~

== octave ==

+ ## full comparison of range hex forms
+ h1 = num2hex (range_test1);
+ h2 = num2hex (range_test2);
+ h3 = num2hex (range_test3);
+ h4 = num2hex (range_test4);
+ sum (sum (h1 != h2))
ans =          0
+ sum (sum (h2 != h3))
ans =      25373
+ sum (sum (h3 != h4))
ans =          5
+ sum (sum (h4 != h1))
ans =      25374
+ sum (sum (h1 != h3))
ans =      25373
+ sum (sum (h2 != h4))
ans =      25374
+ ## full comparison of (range*double) hex forms
+ h1 = num2hex (range_test1 * pi);
+ h2 = num2hex (range_test2 * pi);
+ h3 = num2hex (range_test3 * pi);
+ h4 = num2hex (range_test4 * pi);
+ sum (sum (h1 != h2))
ans =       9518
+ sum (sum (h2 != h3))
ans =      15053
+ sum (sum (h3 != h4))
ans =          5
+ sum (sum (h4 != h1))
ans =      17261
+ sum (sum (h1 != h3))
ans =      17266
+ sum (sum (h2 != h4))
ans =      15053

== matlab ==

% full comparison ofn range hex forms
h1=num2hex(range_test1);
h2=num2hex(range_test2);
h3=num2hex(range_test3);
h4=num2hex(range_test4);
sum(sum(h1~=h2))
ans =
     0
sum(sum(h2~=h3))
ans =
        6752
sum(sum(h3~=h4))
ans =
       10346
sum(sum(h4~=h1))
ans =
        7968
sum(sum(h1~=h3))
ans =
        6752
sum(sum(h2~=h4))
ans =
        7968

% full comparison of (range*double) hex forms
h1=num2hex(range_test1*pi);
h2=num2hex(range_test2*pi);
h3=num2hex(range_test3*pi);
h4=num2hex(range_test4*pi);
sum(sum(h1~=h2))
ans =
     0
sum(sum(h2~=h3))
ans =
        7275
sum(sum(h3~=h4))
ans =
        9932
sum(sum(h4~=h1))
ans =
        8684
sum(sum(h1~=h3))
ans =
        7275
sum(sum(h2~=h4))
ans =
        8684

----





--
View this message in context: 
http://octave.1599824.n4.nabble.com/A-problem-with-range-objects-and-floating-point-numbers-tp4664900p4666760.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.



reply via email to

[Prev in Thread] Current Thread [Next in Thread]