Floating point accurarcy and COLON

I have just been bitten by an odd floating point comparison issue. At the heart of the issue is that I expected to get true for the following test

x = 0:0.1:0.3;
y = 0:0.1:0.4;
isequal(x, y(1:end-1))

Instead

 x - y(1:end-1)
 ans =
   1.0e-16 *
         0         0   -0.2776   -0.5551

The help for COLON says J:D:K is the same as [J, J+D, ..., J+m*D] where m = fix((K-J)/D) which I think supports my expectation. It is not quite right since for 0:0.1:0.3, fix((K-J)/D) is 2. I thought possibly COLON was just doing a linspace, but I also get false for

 z = linspace(0, 0.3, 4);
 isequal(x, z)

Interestingly, for all the cases I have tried I get true for

a = 0;
b = 0.3;
c = 0.1;
x = a:c:b;
isequal(x(end), b)

suggesting COLON somehow fixes the ends, but not in the same way that LINSPACE does.

How does COLON work?

[Walter Roberson]

I think it might be Index = Index+D as you go along. This loses resolution differently than Index = J + N*D.

[Sean de Wolski]

Note, you can see the code for linspace

edit linspace

[Daniel Shub]

@Walter shouldn't methods be true for my first example and not end exactly at "b" in my last example?

@Sean COLON is not using LINSPACE

Answer 1

[Matt Fig]

This explains how the colon operator works. It is an interesting read.

http://www.mathworks.com/support/solutions/en/data/1-4FLI96/index.html?solution=1-4FLI96

Here is a relevant portion:

"To counteract such error accumulation, the algorithm of the COLON operator dictates that:

1. The first half of the output vector (in this example ‘v’) is calculated by adding integer multiples of the step ‘d’ to the left-hand endpoint ‘a’. 2. The second half is calculated by subtracting multiples of the step ‘d’ from the right-hand endpoint."

[Daniel Shub]

I had just found it myself. I did help colon and not doc colon. I think the choice makes sense. Yet another learning experience that floats aren't exact.

[Daniel Shub]

The link no longer works, but is available at http://web.archive.org/web/20120213165003/http://www.mathworks.com/support/solutions/en/data/1-4FLI96/index.html? and http://www.mathworks.co.uk/support/solutions/attachment.html?resid=1-4FLIOY&solution=1-4FLI96.

[Jeffrey]

As Daniel mentioned, the technical solution (1-4FLI96) linked in a previous comment is no longer available at that url.

The solution was republished to MATLAB Central Answers here:

http://www.mathworks.com/matlabcentral/answers/143255-how-does-the-colon-operator-work

[Royi Avital]

Jeffrey, Your link isn't available as well.

[Walter Roberson]

Royi, you are right. I will see if I can get that restored.

[Walter Roberson]

As mentioned before, that technical support solution is gone, and so is the Answers page that republished it.

Fortunately the Technical Support Solution was captured by the WayBack Machine at: https://web.archive.org/web/20121120071523/http://www.mathworks.com/support/solutions/en/data/1-4FLI96/index.html?solution=1-4FLI96

Quoting from there:

Date Last Modified: Friday, April 1, 2011
Solution ID:       1-4FLI96
Product:       MATLAB
Reported in Release:       No Release
Platform:       Windows
Operating System:       Windows XP SP2

Subject:

How does the COLON operator work?

Problem Description:

I want to know how the COLON operator works:

For example, in the command

 v = a:d:b

I would like to know if the vector 'v' is calculated by repeated addition of multiples of 'd' to a.

Solution:

v = a:d:b is not constructed using repeated addition as that might accumulate round-off errors, as described in the related Solution below: "How do I determine if the error in my answer is the result of round-off error or a bug?".

In fact, this kind of round-off error is inherent to software or hardware systems that perform calculations using floating point numbers. Since floating point numbers have to be represented within a limited number of bits, the result of some calculations is rounded off in order to fit the result into the finite representation (sign, exponent, mantissa) for floating-point numbers.

An example of a number that cannot be represented exactly in binary floating point is the number 0.1, which, while finite in decimal represntation, would require an infinite number of bits to be exactly represented in floating-point notation. For instance, executing the following in MATLAB:

 0.1+0.1+0.1 == 0.3

will not return true due to such round-off errors.

To counteract such error accumulation, the algorithm of the COLON operator dictates that:

1. The first half of the output vector (in this example ‘v’) is calculated by adding integer multiples of the step ‘d’ to the left-hand endpoint ‘a’.

2. The second half is calculated by subtracting multiples of the step ‘d’ from the right-hand endpoint.

The COLON operator ensures symmetry of the output vector about its midpoint and also that the round off error is minimized at both endpoints. For example, even though 0.01 is not exactly represented in binary,

 v = -1 : 0.01 : 1

consists of 201 floating points numbers centered symmetrically about zero.

For more detailed information on how the COLON operator works, refer to the attached function (COLONOP). This file mimics closely the behavior of the COLON operator.