MATLAB: Floating point accurarcy and COLON

documentationfloating point

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?

Best Answer

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."

Related Solutions

MATLAB: Floating point arithmetics, mathlab’s interpretation of a . after a number

"Why doesn't e yield 0?"

Because most of your values cannot be exactly represented using binary floating point numbers.

The values stored in your computer's memory are only approximations of those values (to a finite precision), and that approximation means that they include some floating point error. So when you perform arithmetic operations on those values, sometimes that floating point error accumulates during the operation.. .and in some cases, that accumulation can cancel-out completely (important note: this cancelation does not change the precision of the output!). This is exactly what your two examples show: one of the operations happens to accumulate some error which remains in the output, another operation happens to cancel out the error so that the output looks like it is "exact" (in fact in both cases the output precision is much the same).

Lets have a look at the exact values that your computer is calculating with:

>> num2strexact(x)
ans =
0.16666666666666662965923251249478198587894439697265625
>> num2strexact(y)
ans =
9.999999999999997779553950749686919152736663818359375e-2

Learn about binary floating point numbers and how they behave:

https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero

https://www.mathworks.com/matlabcentral/answers/316889-why-two-equal-numbers-are-not-equal

http://www.mathworks.com/matlabcentral/answers/102419-how-do-i-determine-if-the-error-in-my-answer-is-the-result-of-round-off-error-or-a-bug

http://www.mathworks.com/matlabcentral/answers/140656-about-numerical-difficulty-of-equality-comparison-operation-between-two-double-real-numbers-in-matlab

This is worth reading as well:

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

http://faculty.tarleton.edu/agapie/documents/cs_343_arch/papers/1991_Goldberg_FloatingPoint.pdf

MATLAB: Array elements not equal despite only difference being their respective end points

I received an answer from Tech Support which clarifies this issue:

This is a very interesting question. This is an expected behaviour, to understand why this happens we need to understand how the colon operator works.

v = a:d:b is not constructed using repeated addition (as one would think) as that might accumulate round-off errors.

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 representation, 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:

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).

function v = colonop(a,d,b)
% COLONOP  Demonstrate how the built-in a:d:b is constructed.
%




%   v = colonop(a,b) constructs v = a:1:b.
%   v = colonop(a,d,b) constructs v = a:d:b.
%
%   v = a:d:b is not constructed using repeated addition.  If the
%   textual representation of d in the source code cannot be
%   exactly represented in binary floating point, then repeated
%   addition will appear to have accumlated roundoff error.  In
%   some cases, d may be so small that the floating point number
%   nearest a+d is actually a.  Here are two imporant examples.
%
%   v = 1-eps : eps/4 : 1+eps is the nine floating point numbers
%   closest to v = 1 + (-4:1:4)*eps/4.  Since the spacing of the
%   floating point numbers between 1-eps and 1 is eps/2 and the
%   spacing between 1 and 1+eps is eps,
%   v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].
%
%   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.
%
%   Ideally, in exact arithmetic, for b > a and d > 0,
%   v = a:d:b should be the vector of length n+1 generated by
%   v = a + (0:n)*d where n = floor((b-a)/d).
%   In floating point arithmetic, the delicate computatations
%   are the value of n, the value of the right hand end point,
%   c = a+n*d, and symmetry about the mid-point.
if nargin < 3
   b = d;
   d = 1;
end
tol = 2.0*eps*max(abs(a),abs(b));
sig = sign(d);
% Exceptional cases.
if ~isfinite(a) || ~isfinite(d) || ~isfinite(b)
  v = NaN;
  return
elseif d == 0 || a < b && d < 0 || b < a && d > 0
  % Result is empty.
  v = zeros(1,0);
  return
end
% n = number of intervals = length(v) - 1.
if a == floor(a) && d == 1
  % Consecutive integers.
  n = floor(b) - a;
elseif a == floor(a) && d == floor(d)
  % Integers with spacing > 1.
  q = floor(a/d);
  r = a - q*d;
  n = floor((b-r)/d) - q;
else
  % General case.
  n =  round((b-a)/d);
  if sig*(a+n*d - b) > tol
     n = n - 1;
  end
end
% c = right hand end point.
c = a + n*d;
if sig*(c-b) > -tol
  c = b;
end
% v should be symmetric about the mid-point.
v = zeros(1,n+1);
k = 0:floor(n/2);
v(1+k) = a + k*d;
v(n+1-k) = c - k*d;
if mod(n,2) == 0
  v(n/2+1) = (a+c)/2;
end

Best Answer

Related Solutions

MATLAB: Floating point arithmetics, mathlab’s interpretation of a . after a number

MATLAB: Array elements not equal despite only difference being their respective end points

Related Question