Chapter 2 Floating points arithmetic
2.1 What should we care?
This Chapter is only an introduction to floating-points arithmetic. The goal is only to make you realize that most numbers that we use in our numerical projects are only approximations. Our results are therefore subject to rounding errors that may accumulate if we are not careful. If you want a complete and detailed presentation of floating-point numbers, see Goldberg (1991) \href{http://neo.dmcs.p.lodz.pl/ak/IEEE754_article.pdf. As an example, the number 0.1 cannot be represented exactly by computers. We will see why in next section. We can verify that by the following small experiment:
## [1] FALSER tells us that 0.1+0.1+0.1 is not equal to 0.3. However, the difference is small as we can see:
## [1] 5.551115e-17However, when we have to deal with thousands of operations, the error may become more substantial. In the next section, we explain briefly how a computer store numbers and why it cannot represent exactly simple numbers such as 0.1.
2.2 What is a floating-point number
A floating-point number is a real number that can be represented exactly by a computer. First, computers store information in binary format (0’s and 1’s). In order to have software that are portable from one computer or operating system to another, we need a stable way or dealing with floating-point numbers. Different results obtain on different computers is enough to make a program crash. The most common standard used in computers is the IEEE 754 (Institute of Electrical and Electronics Engineers). It is a binary standard which means that numbers are represented in base 2. One of the reasons for using the base 2 is related to the tightness of the relative error that comes from the approximation of real numbers by floating-point numbers. Unless we want to build computers, we don’t need to know more about the advantage of using that standard. In binary format, integers are exactly represented by computers: \[\begin{split} 1 & = 1(2^0) =_2 1\\ 2 & = 1(2^1) + 0(2^0)=_2 10 \\ 3 & = 1(2^1)+1(2^0)=_2 11 \\ 4 & = 1(2^2)+0(2^1) + 0(2^0)=_2 100 \\ \vdots &=_2 \vdots\\ N & = \sum_{i=0}^{p} d_i(2^i) =_2 d_pd_{p-1}\cdots d_0 \end{split} \] where \(=_2\) means “equals in base 2” and \(d_i\)’s are either 0 or 1. Of course, the number of bits restricts the number of integers that can be exactly represented. In general, a real floating-point number, \(x\), is represented as follows: \[ x = \pm \left[d_0 + d_1(2^{-1}) + d_2(2^{-2}) + \cdots + d_{p-1}(2^{-(p-1)})\right] 2^e, \] or simply \[ x = \pm \left[d_0.d_1d_2\cdots d_{p-1}\right]_2 2^e, \] where \([]_2\) means that the inside of the brackets is expressed in base 2, \(p\) is called the precision (or the number of digits) and \(e\) the exponent. The term \(d_0.d_1...d_{p-1}\) is called the significant. For example: \[ 4.5 = 2^2+2^{-1} = (1+2^{-3})2^2 = [1.001]_2\times 2^2, \] which implies that \(d_0=1\), \(d_3=1\) (or a significant equals to 1.001), \(e=2\) (or 10 in base 2) and all other \(d_i\) are equal to zero. Of course there are more than one ways to represent 4.5. Here is another one: \[ 4.5 = (2^{-1}+2^{-4})2^3, \] which implies \(d_1=d_4=1\), \(e=3\) and all other \(d_i\)’s equal to zero. However, the number is uniquely represented if we impose the normalization \(d_0=1\). Furthermore, the normalization saves us a bit because there is no need to store \(d_0\) since it is always 1. The bits are allocated between the storage of the exponent, the \((p-1)\) \(d_i\)’s and the sign (e.g. 1 if + and 0 if -). The exponent is store as an unsigned integer, which implies that the number of different exponents is \((e_{max}-e_{min}\)). Therefore, the number of bits required to store a floating-point number is approximately: \[ Bits = \log_2{(e_{max}-e_{min})}+ (p-1) + 1 \] In R the allocation is stored in the variable .Machine:
## [1] 53## [1] 1024## [1] -1022## [1] 63.99859It gives 1 bit for the sign, 11 for the exponent and 52 for the significant. In fact, R follows the IEEE-754 standard for double precision floating-point numbers. Other important numbers are stored in .Machine. The first is the machine-epsilon (\(\epsilon\)). It is usually defined as \(\beta^{-p+1}/2\), where \(\beta\) is the base. It is the smallest positive floating-point number, \(x\), such that \((1+x)\neq 1\):
## [1] 2.220446e-16## [1] FALSEIt is not the smallest number, which is rather equal to \(2^{e_{min}}\), but the highest relative error that comes from approximating a real numbers by its nearest floating-point number. The accuracy of algorithms is often measured in terms of the machine-epsilon. A bad algorithms may produce errors of size \(n\epsilon\) with \(n>>1\). We therefore want to choose algorithms that minimizes \(n\). In the example we showed at the beginning of the chapter, R returned FALSE to the question: \(0.1+0.1+0.1=0.3\)? But, since 0.1 does not have an exact floating-point representation there is a rounding error that accumulates when summing them. If we want to compare the numbers we need to take into account these rounding errors. The function all.equal compares numbers with a certain level of tolerance. So, they don’t necessarily have to be exactly equal. If we go back to are first example:
## [1] FALSE## [1] TRUEThe interpretation of the above result is that \((.1+.1+.1)\neq 0.3\) only because of rounding errors. The maximum and minimum floating-point numbers are:
## [1] 1.797693e+308## [1] 2.225074e-308Anything above the xmax is considered to be equal to infinity an anything below xmin is considered to be equal to 0:
## [1] InfWe conclude this section by looking at another important results that you may obtain when performing some operations. They are NA or NaN, which are undetermined (or missing) numbers. It will result from operations that produced indeterminate results (or a missing observation when loaded from a dataset). It is not to be confused with operations that produce infinity. For example, \(1/0=\infty\) but \(0/0=NaN\):
## [1] Inf## [1] NaNHere are other examples:
## [1] NaN## Warning in log(-2): NaNs produced## [1] NaN## Warning in sqrt(-1): NaNs produced## [1] NaN## [1] NaNOf course, we cannot compare undetermined results. Therefore, we cannot use the logical operator “==” to verify if an operation produces an NA. In fact the result will also be an NA. If we want to verify if an operation produces an NA or NaN, we use the function is.na:
## [1] NA## [1] TRUEHowever, we can use logical operators with infinity:
## [1] TRUE## [1] TRUE## [1] TRUE2.3 Rounding errors
This section briefly covers the issue of rounding errors created by numerical computations. Again, we cannot cover it in details. In fact, it may be quite complicated to measure rounding errors from certain algorithms. But knowing how they are created may help us developing good programming habits. It will also help justifying some methods that will be presented in the next chapters.
In order to illustrate the problem related to rounding errors, consider the case of floating-point numbers represented in base 10 with the number of digits \(p\) equals to 3 (remember that the number of digits after the point is equal to \((p-1)\)). In that case, \(\pi\) is approximated by the floating-point number \(3.14\times 10^0\). Suppose the true value is 3.1416. The term “last place” refers to the last decimal given by the floating-point representation (or \(10^{-p+1}\)). The error is \(0.0016=0.16\times 10^{-2}\) or \(0.16ulps\) (units of the last place). An error smaller than \(1ulps\) means that the last digit is not contaminated. Error from approximating a real number is always less than \(1ulps\). However, we often get errors greater than \(1ulps\) when results come from mathematical operations. As an example, suppose the result of a computation is \(3.12\times 10^0\) and that the true answer is \(3.14\times 10^0\). The error is then equal to \(2ulps\). In the IEEE standard, where \(p=53\), the precision is up to about the \(16^{th}\) decimal (the \(53^{th}\) digit in base 2 refers to the \(16{th}\) digit in base 10).
The error measure in \(ulps\) is affected by multiplications. Consider the approximated \(\pi\) above. Since it is represented by \(3.14\times 10^0\), if we multiply it by 2, the result will become \(6.28\times 10^0\). The true value is \(2\pi=6.2832\) which implies an error of \(0.32ulps\). Another measure of rounding error is the relative error defined as the error divided by the true value and it is often expressed in term of the machine-epsilon. The relative error for \(\pi\) is \(0.0016/3.1416=0.0005\) and the one for the computation of \(2\pi\) is \(0.0032/6.2832=0.0005\). It is therefore unaffected by the operation. The relative error in terms of the machine-epsilon is \([error/(TRUE*\epsilon)]\epsilon\). In our example, the machine-epsilon is \(10^{-2}/2=0.005\). The relative error for \(\pi\) or \(2\pi\) is therefore \(0.1\epsilon\).
There is a link between the error and the number of contaminated digits. If the error is \(n~ulps\) the number of contaminated digits is \(\log_\beta{n}\), if the relative error is \(n\epsilon\) the number of contaminated digits is \(\log_\beta{n}\). Suppose we want to compute \((x-y)\). Suppose also that the base is 10, p=3 and the computer only keeps 2 decimal when performing operations. In the first case, \(y=2.15\times 10^{12}\) and \(x=1,25\times 10^{-5}\). The computer rewrites the numbers using the same exponent and keep only 2 digits: \[ y-x = 2.15\times 10^{12} - 0.00 \times 10^{12} = 2.15\times 10^{12} \] The error is very small and no digits are contaminated. Consider the second case for which \(y=10.1\) and \(x=9.93\). Using the same rule, we obtain: \[ y-x=1.01\times 10^1-0.99\times 10^1 = 2.00\times 10^{-1} \] The error is \(30 ulps\) which implies a relative error of \([0.03/(0.2*0.005)]\epsilon=30\epsilon\). \(\log_{10}{30}>1\) meaning that the two digits are contaminated. The error from the last subtraction is called “Catastrophic cancellation”, and it arises when \(x\) and \(y\) are very close to each other. Of course, we would not have obtained the same error with double precision floating-point numbers (\(p=53\) and \(\beta=2\)). This example is meant to easily show what happens when we subtract numbers of similar values. In reality, computers are smarter than that. They usually have what we call “Guard Digits” which are extra digits used during floating-point operations. In our example, one Guard Digit would have been enough to have no error. However, “Catastrophic cancellation” exists even in double precision systems. The error is small in absolute term but quite big if we compare it with the error from other floating-point operations. To have an idea in double precision computations, consider the following example:
## [1] 4.421963e-14The error is therefore approximately \(442ulps\) (if we consider that p=53 corresponds to the \(16^{th}\) decimals). However, what we just computed is not really the error coming from \((y-x)\). In fact, rounding errors is hard to measure with a computer because we are not really subtracting the true solution (0.000005 does not have an exact floating-point representation). Catastrophic cancellation occurs in fact when we subtract numbers who are subject to rounding errors. In other words, numbers that come from floating-point operations (ex. \(y^2-x^2\)). In that case, the rounding error is greater than the one from approximating a real number. If the numbers are close, the subtraction eliminates the good digits and we are left with the bad ones. Lets see what happens if \(y\) and \(x\) are not close:
## [1] 0We don’t really care in economics about an error at the \(14^{th}\) decimal. We don’t build bridges. Therefore, we won’t worry when the number of floating-point operations is small. But, these small errors can become huge if the operations are repeated hundreds or thousands of times. For example, when we work with big systems of matrices, solving the system requires many sums, products, subtractions and divisions. If we are not careful about how to minimize rounding errors, we may end up with a wrong solution. Even if we do use some accurate algorithms, there exists iterative procedures to reduce the rounding error further. We will cover that in the next section. Consider only the operation to compute the OLS estimate: \[ \hat{\beta} = (X'X)^{-1}X'Y \] Suppose the sample size is 1000 and the number of coefficients is 10. Then, we need to compute 1000 sums and multiplications for each of the 100 elements of \((X'X)\) and each of the 10 elements of \((X'Y)\). But, that’s not the end. We now have to invert \((X'X)\) (we don’t really do that, as we’ll see in the next chapter) and multiply the inverse by (X’Y). Fortunately, there exists a method that is less sensitive to rounding errors. We don’t even have to compute \(X'X\). We’ll cover that in the next chapter.
The problem with floating-point arithmetic is that a method can be accurate most of the time but very bad in some cases, when our variables take some particular values. And because there is no error message that appears, we think that everything is fine. The purpose of this chapter is only to make sure that you are aware of that. It will also justify many numerical approaches that we will cover throughout the course.
We conclude this chapter with two examples. Suppose we want to compute the sum of the element of a vector. If the elements of the vector come from floating-point operations, the usual sum, \(\sum_{i=1}^n x_i\), could produce large rounding errors. It would be even worse if some successive \(x_i\)’s were close to each other in absolute value with different signs, because we would face multiple catastrophic cancellations. The first solution is to order the vector and to sum the element from the smallest to the largest. Adding a small number to a large number may make us loose the last digits of the small one. That’s why sorting may improve the accuracy of the sum. The second is to use the following algorithm:
Kahan Summation Formula
S=X[1]
C=0
for i = 2 \(\to\) n do
Y=X[i]-C
T=S+Y
C=(T-S)-Y
S=T
end for
Return S
Exercise 2.1 Write the R function mySum for the Kahan Algorithm. Create a \(100,000\times 1\) vector with each element equals to 0.1. Compare the accuracy of your function with the R function sum.
The Kahan algorithm may seem to be an obscure way to improve accuracy. It is in fact quite hard to come up with such algorithms. Another way to minimize rounding errors is to minimize the number of operations when it is necessary. since it reduces the accumulation of errors. Suppose for example that we want to compute the following polynomial (see Judd (1998) page 34): \[ S = \sum_{i=0}^n a_ix^i \] The direct approach requires \(n\) multiplications, \(n\) additions and \((n-1)\) exponentiations for a total of \((3n-1)\) operations. By rewriting the sum as (for \(n=4\)): \[ S = a_0 + x(a_1 + x(a_2 + x(a_3 + xa_4))), \] we still have \(n\) additions and \(n\) multiplications but no exponentiation. The algorithm is called the Horner’s method:
Horner Method
(A is (n+1)\(\times\) vector of \(a_i\))
(x is a scalar)
S=A[n+1]
for i = n \(\to\) 1 do
S=A[i]+S*x
end for
Return S
Exercise 2.2 This is Exercise 8 of Chapter 2 from Judd (1998).
a- Write a function to compute a \(n^{th}\) order polynomial using the Horner’s method.
b- Modify the Horner’s method to compute: \[S = \sum_{i=0}^n\sum_{j=0}^n a_{ij} x^iy^j\]
c- Modify the Horner’s method to compute: \[S = \sum_{i=0}^n\sum_{j=0}^n \sum_{l=0}^n a_{ijl} x^iy^jz^l \]