Numerical Hygiene

Introduction

Numerical expressions are vastly used in the field of data science. They are used to perform mathematical operations on data. However sometimes, the results of these operations are not as expected. This is because of the way computers store numbers and perform operations on them. This blog will help you understand how these operations are performed and how to avoid errors.

Floating Point Numbers

Computers store numbers in binary format. This means that the numbers are stored in the form of 0s and 1s. This is called the binary representation of the number. For example, the number 5 is stored as 101 in binary format. The number 101 is a 3 digit number. This means that the computer uses 3 bits to store the number 5. Similarly, the number 10 is stored as 1010 in binary format.

Often, the number of bits used to store a number is fixed. This means that the computer uses the same number of bits to store all numbers. For example, if the computer uses 8 bits to store a number, then the number 5 will be stored as 00000101. This is called the fixed point representation of the number. For float data type, if the computer uses 16 bits to store a number, then the number 5 will be stored as 0100010100000000 and the number 0.001 will be stored as 0001010000011001.

This is called the floating point representation of the number.

• The first bit is used to store the sign of the number. it is used to store the sign of the number. The sign is 0 for positive numbers and 1 for negative numbers.
• The next 5 bits are used to store the exponent of the number. It is used to shift the decimal point of the number.
• The last 10 bits are used to store the mantissa of the number. It is used to store the fractional part of the number.

import bitstring

n = 5

nb = bitstring.BitArray(float=n, length=32)
print(nb.bin, nb.float)

n = 0.001

nb = bitstring.BitArray(float=n, length=16)

print(nb.bin, nb.float)

01000000101000000000000000000000 5.0
0001010000011001 0.0010004043579101562

Here we can see for only 16 bits, the number 5 is stored as 0100010100000000 and the number 0.001 is stored as 0001010000011001. This means that the computer can store only a limited number of digits after the decimal point. This is called the precision of the number. The precision of the number is determined by the number of bits used to store the number. The more bits used to store the number, the higher the precision of the number. The default bits used to store a number in Python is 64 bits.

Rounding Error

When we perform mathematical operations on numbers, the precision of the number decreases. This is called the rounding error. This happens because the computer uses a limited number of bits to store the result of the operation. This means that the computer cannot store the exact result of the operation. It can only store an approximation of the result. This approximation is called the rounding error.

Sometimes, just the order of operations can cause a rounding error. For example, for extremely small positive number \(\left(x\right)\) and \(\left(x + 1 - 1\right)\) is not equal to \(x\). This is because the computer first adds \(1\) to \(x\) and then subtracts \(1\) from the result. This causes a rounding error. However, if we first subtract \(1\) from \(1\) and then add \(x\) to the result, then the result will be equal to \(x\). This was just a very simple example. However, in some cases, the cause of the rounding error is not so obvious. We will see such cases in this blog.

import numpy as np
from scipy.special import xlog1py, xlogy
import plotly.express as px

import warnings
warnings.filterwarnings("ignore")

import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook_connected"

Softplus

In machine learning, we often use the softplus activation function. It is a smooth approximation to the relu function.

It is defined as:

\[\text{softplus}(x) = \log(1 + e^x)\]

The inverse of the softplus function is the softminus function. It is defined as:

\[\text{softminus}(x) = \log(e^x - 1)\]

def basic_softplus(x):
    y = np.log(1 + np.exp(x))
    # if y = 0, then y = exp(-50)
    y = np.where(y > 0, y, np.exp(-50))
    return y

def relu(x):
    return np.maximum(x, 0)

x = np.linspace(-5, 5, 1000)
y1 = (basic_softplus(x))
y2 = (relu(x))

fig = px.line(title="Softplus function")
fig.add_scatter(x=x, y=y1, name="Softplus")
fig.add_scatter(x=x, y=y2, name="ReLU")
fig.show()

For large negative values of \(x\), the softplus function reaches almost \(0\).

For small positive values of \(x\), the softminus function reaches almost \(-\infty\).

def stable_softplus(x):
    return np.log1p(np.exp(x))

x = np.linspace(-38, -33, 1000)
y1 = np.log(basic_softplus(x))
y2 = np.log(stable_softplus(x))

fig = px.line(title="Log Softplus function")
fig.add_scatter(x=x, y=y1, name="Basic softplus")
fig.add_scatter(x=x, y=y2, name="Numerically stable softplus")
fig.show()

0.3, (3*0.1)

(0.3, 0.30000000000000004)

fig = px.line(title="Difference between Log Softplus and Log Stable Softplus")
fig.add_scatter(x=x, y=(y1 - y2))
fig.show()

Here, we can see that defining the softplus function like basic_softplus causes a rounding error. So, we again defined the softplus function as stable_softplus to avoid the rounding error. We will get more understanding about that in the following sections.

\(\log(1+x)\)

Limit

\[\lim_{x \to 0} \dfrac{x}{\log\left(1+x\right)} = \lim_{x \to 0} \dfrac{1}{\left(\dfrac{1}{1+x}\right)\left(1\right)} = \lim_{x \to 0} \left(1+x\right) = 1\]

Taylor Series Expansion

\[\log(1+x) = \sum_{n=1}^{\infty} \dfrac{(-1)^{n+1}x^n}{n} = x - \dfrac{x^2}{2} + \dfrac{x^3}{3} - \dfrac{x^4}{4} + \dots\]

Therefore, \(\log(1+x) \approx x\) for small \(x\) as the other terms become negligibly tiny.

The log1p function is a useful tool for calculating the logarithm of \(1+x\), particularly when \(x\) is a small value. It helps avoid numerical instability issues that can arise when directly computing np.log(1 + x) when \(x\) is close to zero. By using np.log1p, you can achieve a more accurate and stable calculation of the logarithm.

A = (0, 1, 2, 3, 4)

def logp1(x):
    if np.isinf(x) and x > 0: 
        return x

    u = 1. + x
    d = u - 1.
    
    if d == 0:
        return x
    
    return np.log(u) * x / d

for a in A:
    x = 10**(-13 - a)

    print(f"p = {x}")
    print(f"1 + p = {1 + x}")
    print(f"p - log(1 + p) = {x - np.log(1 + x)}")
    print(f"p - log1p(p) = {x - np.log1p(x)}")
    print(f"p - logp1(p) = {x - logp1(x)}")
    print()

p = 1e-13
1 + p = 1.0000000000001
p - log(1 + p) = 7.992778374091262e-17
p - log1p(p) = 4.998222695480331e-27
p - logp1(p) = 4.998222695480331e-27

p = 1e-14
1 + p = 1.00000000000001
p - log(1 + p) = 7.992778373641611e-18
p - log1p(p) = 5.048709793414476e-29
p - logp1(p) = 5.048709793414476e-29

p = 1e-15
1 + p = 1.000000000000001
p - log(1 + p) = -1.1022302462515587e-16
p - log1p(p) = 5.9164567891575885e-31
p - logp1(p) = 5.9164567891575885e-31

p = 1e-16
1 + p = 1.0
p - log(1 + p) = 1e-16
p - log1p(p) = 0.0
p - logp1(p) = 0.0

p = 1e-17
1 + p = 1.0
p - log(1 + p) = 1e-17
p - log1p(p) = 0.0
p - logp1(p) = 0.0

Here we can see that for a very small value of \(x\), np.log(1 + x) gives result as 0 which is not correct. But, np.log1p(x) gives the almost correct result as x.

\(e^x - 1\)

Limit

\[\lim_{x \to 0} \dfrac{x}{e^x - 1} = \lim_{x \to 0} \dfrac{1}{\left(e^x\right)\left(1\right)} = \lim_{x \to 0} \dfrac{1}{e^x} = 1\]

Taylor Series Expansion

\[e^x = \sum_{n=0}^{\infty} \dfrac{x^n}{n!} = 1 + x + \dfrac{x^2}{2} + \dfrac{x^3}{6} + \dfrac{x^4}{24} + \dots\] \[e^x - 1 = \sum_{n=1}^{\infty} \dfrac{x^n}{n!} = x + \dfrac{x^2}{2} + \dfrac{x^3}{6} + \dfrac{x^4}{24} + \dots\]

Therefore, \(e^x - 1 \approx x\) for very small \(x\) as the other terms become negligibly tiny.

When dealing with small values, using the expression np.exp(x) - 1 can result in a loss of precision. In such cases, the function np.expm1(x) is a better alternative as it provides a more accurate result that preserves the precision of the small input value \(x\).

def exp1m(x):
    if np.isinf(x) and x > 0:
        return x

    u = np.exp(x)
    d = u - 1.

    if d == 0:
        return x
    
    return d * x / np.log(u)

for a in A:
    x = 10**(-13 - a)

    print(f"p = {x}")
    print(f"exp(p) - 1 = {np.exp(x) - 1}")
    print(f"p - (exp(p) - 1) = {x - (np.exp(x) - 1)}")
    print(f"p - expm1(p) = {x - np.expm1(x)}")
    print(f"p - exp1m(p) = {x - exp1m(x)}")
    print()

p = 1e-13
exp(p) - 1 = 9.992007221626409e-14
p - (exp(p) - 1) = 7.99277837359144e-17
p - expm1(p) = -4.998222695480331e-27
p - exp1m(p) = -4.998222695480331e-27

p = 1e-14
exp(p) - 1 = 9.992007221626409e-15
p - (exp(p) - 1) = 7.992778373591124e-18
p - expm1(p) = -5.048709793414476e-29
p - exp1m(p) = -5.048709793414476e-29

p = 1e-15
exp(p) - 1 = 1.1102230246251565e-15
p - (exp(p) - 1) = -1.1022302462515646e-16
p - expm1(p) = -5.9164567891575885e-31
p - exp1m(p) = -5.9164567891575885e-31

p = 1e-16
exp(p) - 1 = 0.0
p - (exp(p) - 1) = 1e-16
p - expm1(p) = 0.0
p - exp1m(p) = 0.0

p = 1e-17
exp(p) - 1 = 0.0
p - (exp(p) - 1) = 1e-17
p - expm1(p) = 0.0
p - exp1m(p) = 0.0

Here we see that for a very small value of \(x\), np.exp(x) - 1 gives result as 0 which is not correct. But, np.expm1(x) gives the almost correct result as x.

\(x \log(y)\)

When \(x = 0\) and we are calculating \(x \log(y)\), we want the result to be \(0\) regardless of the value of \(y\). However, if \(y = 0\), then the result is undefined. Therefore, we need to handle this case separately.

In this case, using the expression x * log(y) gives a result of nan. Here, xlogy(x, y) is a better alternative as it first checks if \(x = 0\) and returns 0 if that is the case. Otherwise, it returns x * log(y).

def xlogyf(x, y):
    print(f"x = {x}, y = {y}")
    print(f"x * log(y) = {x * np.log(y)}")
    print(f"xlogy(y) = {xlogy(x, y)}")
    print()

xlogyf(5, 10)
xlogyf(0, 10)
xlogyf(5, 0)
xlogyf(0, 0)

x = 5, y = 10
x * log(y) = 11.51292546497023
xlogy(y) = 11.51292546497023

x = 0, y = 10
x * log(y) = 0.0
xlogy(y) = 0.0

x = 5, y = 0
x * log(y) = -inf
xlogy(y) = -inf

x = 0, y = 0
x * log(y) = nan
xlogy(y) = 0.0

We can use this xlogy function (and divide by \(\log(2)\)) when calculating the Cross-Entropy Loss.

\(x \log(1+y)\)

Combining two of the above cases,

• When \(x = 0\), we want the result to be \(0\) regardless of the value of \(y\).
• When \(y \approx 0\), we don’t want to lose precision or face numerical instability issues.

In this case, using the expression x * log(1 + y) may result in nan or loss of precision. Here, xlog1py(x, y) is a better alternative for the above reasons.

def xlogyp1(x, y):
    print(f"x = {x}, y = {y}, y + 1 = {y + 1}")
    print(f"x * log(1+y) = {x * np.log(1+y)}")
    print(f"xlog1py(y) = {xlog1py(x, y)}")
    print()

xlogyp1(10, 9)
xlogyp1(1e20, 1e-16)
xlogyp1(0, 1e-16)
xlogyp1(0, 0)
xlogyp1(0, -1)

x = 10, y = 9, y + 1 = 10
x * log(1+y) = 23.02585092994046
xlog1py(y) = 23.02585092994046

x = 1e+20, y = 1e-16, y + 1 = 1.0
x * log(1+y) = 0.0
xlog1py(y) = 10000.0

x = 0, y = 1e-16, y + 1 = 1.0
x * log(1+y) = 0.0
xlog1py(y) = 0.0

x = 0, y = 0, y + 1 = 1
x * log(1+y) = 0.0
xlog1py(y) = 0.0

x = 0, y = -1, y + 1 = 0
x * log(1+y) = nan
xlog1py(y) = 0.0

\(\log(e^{x} + e^{y})\)

The generalisation of \(\log(1+y)\) is \(\log(e^{x} + e^{y})\). Just like the above cases, we want to avoid loss of precision and numerical instability issues. Therefore, instead of using np.log(np.exp(x) + np.exp(y)), we can use np.logaddexp(x, y) as it is a better alternative.

def logpexp(x, y):
    print(f"x = {x}, y = {y}")
    print(f"exp(x) = {np.exp(x)}, exp(y) = {np.exp(y)}")
    print(f"log(exp(x) + exp(y)) = {np.log(np.exp(x) + np.exp(y))}")
    print(f"logaddexp(x, y) = {np.logaddexp(x, y)}")
    print()

logpexp(10, 9)
logpexp(1e20, 1e-16)
logpexp(0, 1e-16)
logpexp(0, 0)
logpexp(0, -1)

x = 10, y = 9
exp(x) = 22026.465794806718, exp(y) = 8103.083927575384
log(exp(x) + exp(y)) = 10.313261687518223
logaddexp(x, y) = 10.313261687518223

x = 1e+20, y = 1e-16
exp(x) = inf, exp(y) = 1.0
log(exp(x) + exp(y)) = inf
logaddexp(x, y) = 1e+20

x = 0, y = 1e-16
exp(x) = 1.0, exp(y) = 1.0
log(exp(x) + exp(y)) = 0.6931471805599453
logaddexp(x, y) = 0.6931471805599454

x = 0, y = 0
exp(x) = 1.0, exp(y) = 1.0
log(exp(x) + exp(y)) = 0.6931471805599453
logaddexp(x, y) = 0.6931471805599453

x = 0, y = -1
exp(x) = 1.0, exp(y) = 0.36787944117144233
log(exp(x) + exp(y)) = 0.31326168751822286
logaddexp(x, y) = 0.31326168751822286