Why Is Colon Operator Not Equivalent To linspace Command?

Why does q not equal zero in the following?

q = isequal(0:.08:50000,linspace(0,50000,625001))

[Voss]

A = 0:0.08:50000;
B = linspace(0,50000,625001);
C = (0:50000*12.5)/12.5;

isequal(A,B)

ans = logical
   0

isequal(A,C)

ans = logical
   0

isequal(B,C)

ans = logical
   1

Answer 1

[James Tursa]

The short answer is that the algorithms are different. In the first place, 0.08 cannot be represented exactly in IEEE double precision floating point, so using that as the delta does result in some binary floating point calculation inaccuracies compared to what you might expect looking at the decimal version. E.g., if you read the doc for colon, it doesn't even guarantee that you get the requested end point if the delta is not an integer. On the other hand, the linspace algorithm always gives you the end points. Also, the intermediate points might not match up exactly either, again because of floating point effects based on the different algorithms. E.g.,

a = 0:.08:50000;
b = linspace(0,50000,625001);
numel(a) == numel(b)

ans = logical
   1

max(abs(a-b))

ans = 7.2760e-12

In your case, you happened to get the same number of elements with both methods, but the points generated are not exactly the same because of floating point effects of the different algorithms.

k = find(a ~= b,1,'last')

k = 573782

eps(a(k))

ans = 7.2760e-12

num2hex(a(k))

ans = '40e669cf5c28f5c2'

num2hex(b(k))

ans = '40e669cf5c28f5c3'

And above you can see the max difference is on the order of eps of the largest element that had a difference. I.e., the very last bit of the floating point value. Generally I would prefer linspace over colon in cases like yours so that you have control over the end points being generated.

Answer 2

[Torsten]

Both arrays seem to be generated differently such that floating point issues come into play:

z1 = 0:.08:50000;
z2 = linspace(0,50000,625001);
[~,idx] = find(z1~=z2);
format long
z1(idx)-z2(idx)

ans = 1×69724

1.0e-11 * 0.000044408920985 0.000044408920985 0.000044408920985 0.000088817841970 0.000088817841970 0.000088817841970 0.000088817841970 0.000088817841970 0.000088817841970 0.000088817841970 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000177635683940 0.000355271367880 0.000355271367880 0.000355271367880 0.000355271367880 0.000355271367880 0.000355271367880

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 3

[Walter Roberson]

The colon operator is defined to use repeated addition. When the increment does not happen to be an integer multiple of a power of 2, there are round-off error problems. For example, 0.08 is not represented by the rational fraction 8 / 100

fprintf('%.999g\n', 0.08)

0.08000000000000000166533453693773481063544750213623046875

those ....1665<etc> accumulate under repeated addition.

The colon operator is not defined as (0:NumberOfValues-1)*Increment+BaseNumber . If it were defined that way, then

for K = 1:inf

would have to generate a vector 1, 2, 3, ... 9223372036854775806 ahead of time and use the elements from that vector. But 9223372036854775806 is 2^63 and it is not possible to generate vectors longer than 2^47-1 even if memory was available.

Of course, potentially MATLAB could have defined that the colon operator in

 K=1:1.1:2^40

worked differently than the colon operator in

 for K=1:1.1:2^40

but it chose to define them the same way: using repeated addition both ways.