Testes numéricos e estabilidade

Igualdade com tolerância

Dois números calculados por caminhos diferentes raramente são idênticos até o último bit. Em vez de ==, compare com tolerância. A regra do isclose é |a − b| <= atol + rtol × |b|, com rtol=1e-5 e atol=1e-8 por padrão. A parte relativa serve para números grandes, e a absoluta, para números perto de 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

Comparando com zero, a tolerância relativa não ajuda (qualquer coisa é infinitamente maior que zero em termos relativos), e é a parte absoluta que decide. Por isso um teste que compara com zero precisa de um atol escolhido conscientemente.

`numpy.testing`

Para testes automatizados, o módulo numpy.testing dá mensagens úteis: em vez de um simples False, ele mostra quais elementos diferem e por quanto:

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

Dentro do pytest, prefira assert_allclose a assert np.allclose(...): quando o teste falha, você vê o motivo, e não apenas "assertion failed".

A ordem da soma muda o resultado

A adição de decimais não é associativa: (a + b) + c pode diferir de a + (b + c). Um número pequeno somado a um enorme "some" (perde-se na precisão), e se o enorme for cancelado depois, o pequeno já foi embora:

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

A soma de esquerda para direita deu 0, e somando primeiro os dois grandes deu 1, o resultado correto, que o math.fsum (uma soma de alta precisão) também acha. O NumPy usa uma soma "em pares" (pairwise) para arrays grandes, que reduz esse erro em comparação com uma soma sequencial:

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

1.6777216e+07 2e+07

O cumsum soma sequencialmente em float32 e trava em 16.777.216 (a partir daí, somar 1 não muda nada, como vimos no capítulo 7). O sum, com a soma em pares, chega ao valor certo de 20 milhões (o Python imprime 2e+07). Conhecer isso evita acusar o NumPy de "errar" uma conta que só é imprecisa no tipo que você escolheu.

Cancelamento catastrófico

Subtrair dois números quase iguais destrói dígitos significativos. Quando existe uma função feita para o caso, use-a. O np.log1p(x) calcula log(1 + x) com precisão mesmo para x minúsculo, e o expm1 faz o mesmo para exp(x) − 1:

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

1.000000082690371e-10 9.999999999500001e-11

O primeiro resultado perdeu quase metade dos dígitos: o 1 + y foi arredondado antes do logaritmo.

Estabilidade numérica: o softmax

Calcular o softmax ingenuamente estoura com valores grandes, porque o exp(1000) é infinito, e inf / inf é nan. Subtrair o máximo antes (que não muda o resultado matemático) resolve:

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]

Verificar um gradiente por diferenças finitas

Quando você escreve o gradiente de uma função à mão (o caso das redes neurais do capítulo 43), é fácil errar um sinal. Uma checagem confiável compara o gradiente analítico com a diferença finita: (f(w + h) − f(w − h)) / 2h, uma estimativa numérica da derivada:

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