Scalars

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

Integers

Binary representation of integers

format(16, '032b')

'00000000000000000000000000010000'

Bit shifting

format(16 >> 2, '032b')

'00000000000000000000000000000100'

16 >> 2

4

format(16 << 2, '032b')

'00000000000000000000000001000000'

16 << 2

64

Endianess

This refers to how the bytes that make up an integer are stored in the computer. Big Endian means that the most significant byte (8 bits) is stored at the lowest address, while Little Endian means that the least significant byte is stored at the lowest address.

Source: https://chortle.ccsu.edu/AssemblyTutorial/Chapter-15/bigLittleEndian.gif

Overflow

In general, the computer representation of integers has a limited range, and may overflow. The range depends on whether the integer is signed or unsigned.

For example, with 8 bits, we can represent at most \(2^8 = 256\) integers.

• 0 to 255 unsigned

• -128 ti 127 signed

Signed integers

np.arange(130, dtype=np.int8)[-5:]

array([ 125,  126,  127, -128, -127], dtype=int8)

Unsigned integers

np.arange(130, dtype=np.uint8)[-5:]

array([125, 126, 127, 128, 129], dtype=uint8)

np.arange(260, dtype=np.uint8)[-5:]

array([255,   0,   1,   2,   3], dtype=uint8)

Integer division

In Python 2 or other languages such as C/C++, be very careful when dividing as the division operator / performs integer division when both numerator and denominator are integers. This is rarely what you want. In Python 3 the / always performs floating point division, and you use // for integer division, removing a common source of bugs in numerical calculations.

%%python2

import numpy as np

x = np.arange(10)
print(x/10)

Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
ImportError: No module named numpy

Python 3 does the “right” thing.

x = np.arange(10)
x/10

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

Real numbers

Real numbers are represented as floating point numbers. A floating point number is stored in 3 pieces (sign bit, exponent, mantissa) so that every float is represetned as get +/- mantissa ^ exponent. Because of this, the interval between consecutive numbers is smallest (high precison) for numebrs close to 0 and largest for numbers close to the lower and upper bounds.

Because exponents have to be singed to represent both small and large numbers, but it is more convenint to use unsigned numbers here, the exponnent has an offset (also knwnn as the exponentn bias). For example, if the expoennt is an unsigned 8-bit number, it can rerpesent the range (0, 255). By using an offset of 128, it will now represent the range (-127, 128).

Note: Intervals between consecutive floating point numbers are not constant. In particular, the precision for small numbers is much larger than for large numbers. In fact, approximately half of all floating point numbers lie between -1 and 1 when using the double type in C/C++ (also the default for numpy).

Because of this, if you are adding many numbers, it is more accurate to first add the small numbers before the large numbers.

See Wikipedia for how this binary number is evaluated to 0.15625.

from ctypes import c_int, c_float

s = c_int.from_buffer(c_float(0.15625)).value

s = format(s, '032b')
s

'00111110001000000000000000000000'

rep = {
    'sign': s[:1],
    'exponent' : s[1:9:],
    'fraction' : s[9:]
}
rep

{'exponent': '01111100', 'fraction': '01000000000000000000000', 'sign': '0'}

Most base 10 real numbers are approximations

This is simply because numbers are stored in finite-precision binary format.

'%.20f' % (0.1 * 0.1 * 100)

'1.00000000000000022204'

Never check for equality of floating point numbers

i = 0
loops = 0
while i != 1:
    i += 0.1 * 0.1
    loops += 1
    if loops == 1000000:
        break
i

10000.000000171856

i = 0
loops = 0
while np.abs(1 - i) > 1e-6:
    i += 0.1 * 0.1
    loops += 1
    if loops == 1000000:
        break
i

1.0000000000000007

Associative law does not necessarily hold

6.022e23 - 6.022e23 + 1

1.0

1 + 6.022e23 - 6.022e23

0.0

Distributive law does not hold

a = np.exp(1)
b = np.pi
c = np.sin(1)

a*(b+c)

10.82708950985241

a*b + a*c

10.827089509852408

Catastrophic cancellation

Consider calculating sample variance

\[s^2= \frac{1}{n(n-1)}\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2\]

Be careful whenever you calculate the difference of potentially big numbers.

def var(x):
    """Returns variance of sample data using sum of squares formula."""

    n = len(x)
    return (1.0/(n*(n-1))*(n*np.sum(x**2) - (np.sum(x))**2))

Numerically stable algorithms

What is the sample variance for numbers from a normal distribution with variance 1?

np.random.seed(15)
x_ = np.random.normal(0, 1, int(1e6))
x = 1e12 + x_
var(x)

-1328166901.4739892

np.var(x)

1.001345504504934

Underflow

np.warnings.filterwarnings('ignore')

np.random.seed(4)
xs = np.random.random(1000)
ys = np.random.random(1000)
np.prod(xs)/np.prod(ys)

nan

Prevent underflow by staying in log space

x = np.sum(np.log(xs))
y = np.sum(np.log(ys))
np.exp(x - y)

696432868222.2549

Overflow

Let’s calculate

\[\log(e^{1000} + e^{1000})\]

Using basic algebra, we get the solution \(\log(2) + 1000\).

x = np.array([1000, 1000])
np.log(np.sum(np.exp(x)))

inf

np.logaddexp(*x)

1000.6931471805599

logsumexp

This function generalizes logaddexp to an arbitrary number of addends and is useful in a variety of statistical contexts.

Suppose we need to calculate a probability distribution \(\pi\) parameterized by a vector \(x\)

\[\pi_i = \frac{e^{x_i}}{\sum_{j=1}^n e^{x_j}}\]

Taking logs, we get

\[\log(\pi_i) = x_i - \log{\sum_{j=1}^n e^{x_j}}\]

x = 1e6*np.random.random(100)

np.log(np.sum(np.exp(x)))

inf

from scipy.special import logsumexp

logsumexp(x)

989279.6334854055

Other useful numerically stable functions

logp1 and expm1

np.exp(np.log(1 + 1e-6)) - 1

9.999999999177334e-07

np.expm1(np.log1p(1e-6))

1e-06

sinc

x = 1

np.sin(x)/x

0.8414709848078965

np.sinc(x)

3.8981718325193755e-17

x = np.linspace(0.01, 2*np.pi, 100)

plt.plot(x, np.sinc(x), label='Library function')
plt.plot(x, np.sin(x)/x, label='DIY function')
plt.legend()
pass