Numerical tests and stability

Equality with tolerance

Two numbers calculated by different paths are rarely identical up to the last bit. Instead of ==, compare with tolerance. The rule of isclose is |a − b| <= atol + rtol × |b|, with rtol=1e-5 and atol=1e-8 by default. The relative part serves for large numbers, and the absolute part, for numbers close to zero:

import math

import numpy as np

print(0.1 + 0.2 == 0.3, np.isclose(0.1 + 0.2, 0.3), np.allclose([1.0, 2.0], [1.0 + 1e-9, 2.0]))
print(np.isclose(1e-10, 0.0), np.isclose(1e-10, 0.0, atol=0))

False True True
True False

Comparing with zero, relative tolerance does not help (anything is infinitely larger than zero in relative terms), and it is the absolute part that decides. Therefore, a test that compares with zero needs a atol consciously chosen.

`numpy.testing`

For automated tests, the module numpy.testing provides useful messages: instead of a simple False, it shows which elements differ and by how much:

from numpy.testing import assert_allclose

assert_allclose(np.array([1.0, 2.0]), np.array([1.0 + 1e-9, 2.0]), rtol=1e-6)
try:
    assert_allclose(np.array([1.0, 2.0]), np.array([1.0, 2.1]))
except AssertionError as erro:
    print("Mismatched elements: 1 / 2 (50%)" in str(erro))

True

Within pytest, prefer assert_allclose to assert np.allclose(...): when the test fails, you see the reason, not just "assertion failed".

The order of summation changes the result

The addition of decimals is not associative: (a + b) + c can differ from a + (b + c). A small number added to a huge one "sums" (loses precision), and if the huge one is canceled later, the small one is already gone:

x = np.array([1e16, 1.0, -1e16])
print(x.sum(), x[0] + x[2] + x[1], math.fsum(x))

0.0 1.0 1.0

The sum from left to right gave 0, and summing the two large ones first gave 1, the correct result, which the math.fsum (a high-precision sum) also finds. NumPy uses a "pairwise" sum (pairwise) for large arrays, which reduces this error compared to a sequential sum:

grande = np.ones(20_000_000, dtype=np.float32)
print(np.cumsum(grande)[-1], grande.sum())

1.6777216e+07 2e+07

The cumsum sums sequentially in float32 and locks at 16,777,216 (from there, adding 1 changes nothing, as we saw in chapter 7). The sum, with the sum in pairs, reaches the correct value of 20 million (Python prints 2e+07). Knowing this prevents blaming NumPy for "making a mistake" in a calculation that is only imprecise in the type you chose.

Catastrophic cancellation

Subtracting two nearly equal numbers destroys significant digits. When there is a function made for the case, use it. The np.log1p(x) calculates log(1 + x) accurately even for x tiny, and the expm1 does the same for exp(x) − 1:

y = 1e-10
print(np.log(1 + y), np.log1p(y))

1.000000082690371e-10 9.999999999500001e-11

The first result lost almost half of the digits: the 1 + y was rounded before the logarithm.

Numerical stability: the softmax

Calculating softmax naively overflows with large values because the exp(1000) is infinite, and inf / inf is nan. Subtracting the maximum beforehand (which does not change the mathematical result) solves it:

z = np.array([1000.0, 1001.0, 1002.0])
with np.errstate(over="ignore", invalid="ignore"):
    ingenuo = np.exp(z) / np.exp(z).sum()
estavel = np.exp(z - z.max()) / np.exp(z - z.max()).sum()
print(ingenuo, estavel.round(4))

[nan nan nan] [0.09   0.2447 0.6652]

Verifying a gradient by finite differences

When you write the gradient of a function by hand (the case of neural networks in chapter 43), it's easy to get a sign wrong. A reliable check compares the analytical gradient with the finite difference: (f(w + h) − f(w − h)) / 2h, a numerical estimate of the derivative:

def f(w):
    return (w ** 3).sum() + 2 * w[0] * w[1]


def gradiente_analitico(w):
    g = 3 * w ** 2
    g[0] += 2 * w[1]
    g[1] += 2 * w[0]
    return g


def gradiente_numerico(funcao, w, h=1e-6):
    g = np.zeros_like(w)
    for i in range(w.size):
        mais = w.copy()
        mais[i] += h
        menos = w.copy()
        menos[i] -= h
        g[i] = (funcao(mais) - funcao(menos)) / (2 * h)
    return g


w = np.array([1.0, -2.0, 0.5])
print(np.allclose(gradiente_analitico(w), gradiente_numerico(f, w), atol=1e-5))

True