In [1]:

#
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib as mp
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn
import laUtilities as ut
import slideUtilities as sl
import demoUtilities as dm
from matplotlib import animation
from importlib import reload
from datetime import datetime
from IPython.display import Image, display_html, display, Math, HTML;
qr_setting = None

mp.rcParams['animation.html'] = 'jshtml';

Announcements¶

Test 1 is in-class on Friday
Upcoming office hours:
- Today: Prof McDonald from 4:30-6pm in CCDS 1341
- Tomorrow: Peer tutor Rohan Anand from 12-3pm in CCDS 16th floor
Reading:
- Boyd-Vandenberghe Chapter 6
- 3Blue1Brown video 3 and video 4

Lecture 15: Computing with Real Numbers¶

One of the themes of this course will be shifting between mathematical and computational views of various concepts.

Today we need to talk about why the answers we get from computers can be different from the answers we get mathematically -- for the same question!

The root of the problem has to do with how numbers are manipulated in a computer.

In a computer we assign bit patterns to correspond to certain numbers.

We say the bit pattern is the number's representation.

For example the number '3.14' might have the representation '01000000010010001111010111000011'.

For reasons of efficiency, we use a fixed number of bits for these representations. In most computers nowadays we use 64 bits to represent a number.

15.1 Floating-Point Representations of Real Numbers¶

Representing a real number in a computer is very complicated, much more so than representing integers (where computer arithmetic is, mostly, correct). For many decades after electronic computers were developed, there was no accepted "best" way to do this!

Eventually (in the 1980s) a widely accepted standard emerged, called IEEE-754. This is what almost all computers now use.

The style of representation used is called floating point.

Conceptually, it is similar to "scientific notation."

123456 = \underset{significand}{\underset{⏟}{1.23456}} \times {\underset{base}{\underset{⏟}{10}}}^{\overset{e x p o n e n t}{\overset{⏞}{5}}}

Except that it is encoded in binary:

17 = \underset{s i g n i f i c a n d}{\underset{⏟}{1.0001}} \times {\underset{b a s e}{\underset{⏟}{2}}}^{\overset{e x p o n e n t}{\overset{⏞}{4}}}

In a double-precision floating point number: the sign, significand, and exponent are all contained within 64 bits (or 8 bytes).

In [54]:

#
display(Image("images/04-floating-point.png", width=600))

By Codekaizen (Own work) [GFDL or CC BY-SA 4.0-3.0-2.5-2.0-1.0], via Wikimedia Commons

Special Values¶

There are three kinds of special values defined by IEEE-754:

NaN, which means "Not a Number"
Infinity -- both positive and negative
Zero -- both positive and negative.

NaN and Inf behave about as you'd expect. If you get one of these values in a computation you should be able to reason about how it happened. Note that these are values, and can be assigned to variables.

In [41]:

np.sqrt(-1)

<ipython-input-41-597592b72a04>:1: RuntimeWarning: invalid value encountered in sqrt
  np.sqrt(-1)

Out[41]:

nan

In [42]:

var = np.log(0)
var

<ipython-input-42-49e412a30c50>:1: RuntimeWarning: divide by zero encountered in log
  var = np.log(0)

Out[42]:

-inf

In [8]:

1/var

Out[8]:

-0.0

As far as we are concerned, there is no difference between positive and negative zero. You can ignore the minus sign in front of a negative zero.

In [9]:

var = np.nan
var + 7

Out[9]:

nan

In [10]:

var = np.inf
var + 7

Out[10]:

inf

The Relative Error of a Real Number stored in a Computer¶

Because only a fixed number of bits are used, most real numbers cannot be represented exactly in a computer.

Another way of saying this is that, usually, a floating point number is an approximation of some particular real number.

Floating point numbers cannot be arbitrarily large or small. Numbers can be as large as $1.79 \times 10^{308}$ and as small as $2.23 \times 10^{- 308}$ .
There are discrete gaps between floating point numbers. For instance, the interval $[1, 2]$ is represented by discrete subset: $1, 1 + 2^{- 52}, 1 + 2 \times 2^{- 52}, 1 + 3 \times 2^{- 52}, \dots, 2$

Floats are discrete but not equidistant:

In [52]:

#
display(Image("images/04-scale.png", width=600))

Source: What you never wanted to know about floating point but will be forced to find out

Generally when we try to store a real number in a computer, what we wind up storing is the closest floating point number that the computer can represent.

The way to think about working with floating point (in fact, how the hardware actually does it) is:

Represent each input as the nearest representable floating point number.
Compute the result exactly from the floating point representations.
Return the nearest representable floating point number to the result.

What does "nearest" mean? Long story short, it means "round to the nearest representable value."

Let's say we have a particular real number $r$ and we represent it as a floating point value $f$ .

Then $r = f + ϵ$ where $ϵ$ is the amount that $r$ was rounded when represented as $f$ .

How big can $ϵ$ be? Let's say $r$ is a number that is ever so slightly larger than 1 (or more generally, any power of two $2^{n}$ ). Our computer has two choices: round it to 1 or the next larger number.

f = \underset{53 b i t s}{\underset{⏟}{1.000...00}} \times 2^{n} or f + ϵ = \underset{53 b i t s}{\underset{⏟}{1.000...01}} \times 2^{n}

Then $| ϵ |$ must be smaller than

| ϵ | < \underset{53 b i t s}{\underset{⏟}{0.000...01}} \times 2^{n}

So as a relative error,

relative error = \frac{| ϵ |}{f} < \frac{0.000...01 \times 2^{n}}{\underset{53 b i t s}{\underset{⏟}{1.000...00}} \times 2^{n}} = 2^{- 52} \approx 10^{- 16}

Numpy will tell you the epsilon for your computer. It corresponds to the distance between 1 and the next larger number.

In [35]:

print(np.finfo(float).eps)
print(np.finfo(float).eps == 2**-52)

2.220446049250313e-16
True

This value $10^{- 16}$ is an important one to remember.

It is approximately the relative error that can be introduced any time a real number is stored in a computer.

Another way of thinking about this is that you only have about 16 digits of accuracy in a floating point number.

One implication is that if you add a large value with a small value, the small one might disappear entirely. This is called an underflow error.

In [1]:

a = float(10**10)
b = float(10**-10)
print(a)
print(b)
print(a - b)
a - b == a

10000000000.0
1e-10
10000000000.0

Out[1]:

True

Another implication is that every operation introduces error. If we let $*$ denote any mathematical operation (e.g., addition, subtraction, multiplication, or division) and $⊛$ denote its floating point analogue, $x ⊛ y = (x * y) (1 + ϵ)$ for some $ϵ$ , where $| ϵ | \leq ϵ$ . That is, every operation of floating point arithmetic is exact up to a relative error of size at most $ε$ .

As a result, it's possible to create calculations that give very wrong answers, particularly when repeating an operation many times, since each operation could amplify the error. We saw an example of this a few lectures ago.

In [24]:

def f(x):
    if x <= 1/2:
        return 2 * x
    if x > 1/2:
        return 2*x - 1

x = 0.2
v = np.zeros(80)
for i in range(80):
    v[i] = x
    x = f(x)

print("start:", v[:4])
print("end:  ", v[-4:])

start: [0.2 0.4 0.8 0.6]
end:   [1. 1. 1. 1.]

15.2 Principles for Approximate Computing¶

What does all of this mean for us in practice? Since we cannot represent numbers exactly on a computer, it becomes important to know how small perturbations in the input to a problem impact the output.

Here are 3 principles to keep in mind when computing with floating point numbers.

Principle 1: Compare floating point numbers for closeness, not equality¶

Two floating point computations that should yield the same result mathematically, may not do so due to rounding error.

However, in general, if two numbers should be equal, the relative error of the difference in the floating point should be small.

So, instead of asking whether two floating numbers are equal, we should ask if they are close -- meaning that the relative error of their difference is small. Let's try an example.

In [31]:

a = 7
b = 1/10
c = 1/7
r1 = a * b * c
r2 = b * c * a
print(r1)
print(r2)
np.abs(r1-r2)/r1

0.1
0.09999999999999999

Out[31]:

1.3877787807814457e-16

In [32]:

# The wrong way to check equality
print(r1 == r2)

False

In this case, both 1/70 and 0.1 cannot be stored exactly in a floating point representation.

More importantly, the rounding errors are different for the two numbers. As a result, the process of rounding 1/70 to its closest floating point representation, then multiplying by 7, yields a number whose closest floating point representation is not 0.1

However, that floating point representation is very close to 0.1.

Let's look at the difference: -1.3877787807814457e-17.

This is about $- 1 \cdot 10^{- 17}$ .

In other words, -0.0000000000000001

Compared to 0.1, this is a very small number. The relative error abs(-0.0000000000000001 / 0.1) is about $10^{- 16} .$

This suggests that when a floating point calculation is not exact, the error (in a relative sense) is usually very small.

Notice also that in our example the size of the relative error is about $10^{- 16}$ .

Recall that the significand in IEEE-754 uses 52 bits.

Now, note that $2^{- 52} \approx 10^{- 16}$ .

There's our "sixteen digits of accuracy" principle again.

In [33]:

# A better way to check equality
np.finfo('float')
print(np.abs(r1 -  r2)/np.max([r1, r2]) < np.finfo('float').eps)

True

This test is needed often enough that numpy has a function that implements it:

In [10]:

np.isclose(r1, r2)

Out[10]:

True

The matrix version of this function is called allclose. It applies an isclose check for each entry in the matrix.

In [3]:

A = np.random.rand(2,2)
B = np.random.rand(2,2)
C = np.random.rand(2,2)

(A @ B) @ C == A @ (B @ C)

Out[3]:

array([[False, False],
       [False, False]])

In [4]:

np.allclose((A @ B) @ C, A @ (B @ C))

Out[4]:

True

So another way to state Rule 1 for our purposes is:

… always use np.isclose() or np.allclose() to compare floating point numbers for equality!

Next, we will generalize this idea a bit:

beyond the fact that numbers that should be equal, may not be in practice,

we can also observe that it can be hard to be accurate about the difference between two numbers that are nearly equal. This leads to the next two principles.

Principle 2: Beware of ill-conditioned problems¶

An ill-conditioned problem is one in which the outputs depend in a very sensitive manner on the inputs.

That is, a small change in the inputs can yield a very large change in the outputs.

The simplest example is computing $1 / (a - b)$ . If $a$ is close to $b$ , then small changes in either make a big difference in the output. This is because 0 doesn't have a multiplicative inverse, so as $a$ and $b$ get closer then $1 / (a - b)$ will go to infinity.

In [11]:

print(f'r1 is {r1}')
print(f'r2 is very close to r1')
r3 = r1 + 0.0001
print(f'r3 is 0.1001')
print(f'1/(r1 - r2) = {1/(r1 - r2)}')
print(f'1/(r3 - r2) = {1/(r3 - r2)}')

r1 is 0.1
r2 is very close to r1
r3 is 0.1001
1/(r1 - r2) = 7.205759403792794e+16
1/(r3 - r2) = 9999.999999998327

Because the inputs to your problem may not be exact, if the problem is ill-conditioned, the outputs may be wrong by a large amount.

The notion of ill-conditioning applies to matrix problems too, and in particular comes up when we solve certain problems involving matrices that are "almost uninvertible."

Here is a simple example; we will see later why this issue arises. Given a matrix $A$ , let's try to solve the two equations

A x = [\begin{matrix} 1 \\ 2 \end{matrix}] and A y = [\begin{matrix} 1 \\ 2.01 \end{matrix}] .

The two vectors on the right hand sides are similar, so we might expect that the (one and only) solution $x$ would be close to $y$ . But sometimes this is not the case.

In [48]:

A = np.array([[1, 2.0000000001], [2, 4]])
b = np.array([[1],[2]])
c = np.array([[1],[2.01]])

print("x =")
print(np.linalg.solve(A, b))

print("y =")
print(np.linalg.solve(A, c)) # should be [[100 million], [-50 million]]

x =
[[1.]
 [0.]]
y =
[[ 99999992.73096144]
 [-49999995.86298072]]

Principle 3: Relative error can be magnified during subtractions¶

Two numbers, each with small relative error, can yield a value with large relative error if subtracted. If two numbers have errors in the 10th significant digit, but the first 8 significant digits cancel out, then all of a sudden the relative error is large.

This is easiest to see with an example.

In [16]:

a = 1.23456789
b = 1.2345678
trueDiff = 0.00000009
calcDiff = a - b
print("True difference:", trueDiff)
print("Calculated diff:", calcDiff,"\n")
print("Absolute error:", trueDiff-calcDiff)
print("Relative error:", (trueDiff-calcDiff) / trueDiff)

True difference: 9e-08
Calculated diff: 8.999999989711682e-08 

Absolute error: 1.0288317610723888e-16
Relative error: 1.1431464011915431e-09

Here, the relative error in the inputs is on the order of $10^{- 16}$ , but the relative error of the output is on the order of $10^{- 9}$ – i.e., a million times larger.

A good summary that covers additional issues is at https://docs.python.org/2/tutorial/floatingpoint.html.