Phase 01 - Lesson 13

Estabilidad Numérica

This lesson includes a graded coding exercise that runs in your browser, unlocked with lifetime access.

El punto flotante es una abstracción con fugas. Te va a morder durante el entrenamiento, y no lo vas a ver venir.

Tipo: Construir Lenguaje: Python Requisitos previos: Fase 1, Lecciones 01-04 Tiempo: ~120 minutos

Objetivos de Aprendizaje

• Implementar softmax y log-sum-exp numéricamente estables usando el truco de la resta del máximo
• Identificar overflow, underflow y cancelación catastrófica en computaciones de punto flotante
• Verificar gradientes analíticos contra gradientes numéricos usando diferencias finitas centradas
• Explicar por que bfloat16 se prefiere sobre float16 para el entrenamiento y como el escalado de la pérdida evita el underflow de gradientes

El Problema

Tu modelo entrena durante tres horas y luego la pérdida se vuelve NaN. Agregas una instrucción de print. Los logits están bien en el paso 9.000. En el paso 9.001 son inf. Para el paso 9.002 cada gradiente es nan y el entrenamiento está muerto.

O: tu modelo entrena hasta el final, pero la exactitud es 2% peor de lo que afirma el paper. Revisas todo. La arquitectura coincide. Los hiperparámetros coinciden. Los datos coinciden. El problema es que el paper uso float32 y tu usaste float16 sin el escalado correcto. Treinta y dos bits de error de redondeo acumulado se comieron silenciosamente tu exactitud.

O: implementas la pérdida de entropía cruzada desde cero. Funciona con logits pequeños. Cuando los logits superan 100, devuelve inf. El softmax desbordo porque exp(100) es más grande de lo que float32 puede representar. Todo framework de ML maneja esto con un truco de dos líneas. No sabias que el truco existia.

La estabilidad numérica no es una preocupación teórica. Es la diferencia entre una corrida de entrenamiento que tiene éxito y una que falla silenciosamente. Cada bug serio de ML que llegues a depurar termina reduciéndose a punto flotante.

El Concepto

IEEE 754: Como las Computadoras Almacenan Números Reales

Las computadoras almacenan números reales como valores de punto flotante siguiendo el estándar IEEE 754. Un float tiene tres partes: un bit de signo, un exponente y una mantisa (significando).

Float32 layout (32 bits total):
[1 sign] [8 exponent] [23 mantissa]

Value = (-1)^sign * 2^(exponent - 127) * 1.mantissa

La mantisa determina la precisión (cuantos dígitos significativos). El exponente determina el rango (que tan grande o pequeño puede ser un número).

Format     Bits   Exponent  Mantissa  Decimal digits  Range (approx)
float64    64     11        52        ~15-16          +/- 1.8e308
float32    32     8         23        ~7-8            +/- 3.4e38
float16    16     5         10        ~3-4            +/- 65,504
bfloat16   16     8         7         ~2-3            +/- 3.4e38

float32 te da alrededor de 7 dígitos decimales de precisión. Eso significa que puede distinguir 1.0000001 de 1.0000002, pero no 1.00000001 de 1.00000002. Después de 7 dígitos, todo es ruido de redondeo.

float16 te da alrededor de 3 dígitos. El número más grande que puede representar es 65.504. Eso es inquietantemente pequeño para ML, donde logits, gradientes y activaciones superan ese valor rutinariamente.

bfloat16 es la respuesta de Google al problema de rango de float16. Tiene el mismo exponente de 8 bits que float32 (mismo rango, hasta 3.4e38) pero solo 7 bits de mantisa (menos precisión que float16). Para entrenar redes neuronales, el rango importa más que la precisión, así que bfloat16 normalmente gana.

Por que 0.1 + 0.2 != 0.3

El número 0.1 no se puede representar exactamente en punto flotante binario. En base 2, es una fracción periódica:

0.1 in binary = 0.0001100110011001100110011... (repeating forever)

float32 trunca esto a 23 bits de mantisa. El valor almacenado es aproximadamente 0.100000001490116. De forma similar, 0.2 se almacena como aproximadamente 0.200000002980232. Su suma es 0.300000004470348, no 0.3.

In Python:
>>> 0.1 + 0.2
0.30000000000000004

>>> 0.1 + 0.2 == 0.3
False

Esto importa para ML porque:

1. Las comparaciones de pérdida como if loss < threshold pueden dar respuestas erradas
2. Acumular muchos valores pequeños (actualizaciones de gradiente a lo largo de miles de pasos) se desvía de la suma verdadera
3. Los checksums y las pruebas de reproducibilidad fallan si comparas floats con ==

La corrección: nunca compares floats con ==. Usa abs(a - b) < epsilon o math.isclose().

Cancelación Catastrófica

Cuando restas dos números de punto flotante casi iguales, los dígitos significativos se cancelan y te quedas con ruido de redondeo promovido a dígitos principales.

a = 1.0000001    (stored as 1.00000011920929 in float32)
b = 1.0000000    (stored as 1.00000000000000 in float32)

True difference:  0.0000001
Computed:         0.00000011920929

Relative error: 19.2%

Eso es un error relativo del 19% a partir de una sola resta. En ML, esto pasa cada vez que:

• Calculas la varianza de datos con una media grande: E[x^2] - E[x]^2 cuando E[x] es grande
• Restas log-probabilidades casi iguales
• Calculas gradientes por diferencias finitas con un epsilon demasiado pequeño

La corrección: reorganiza las fórmulas para evitar restar números grandes y casi iguales. Para la varianza, usa el algoritmo de Welford o centra los datos primero. Para las log-probabilidades, trabaja en el espacio logarítmico todo el tiempo.

Overflow y Underflow

El overflow pasa cuando un resultado es demasiado grande para representarse. El underflow pasa cuando es demasiado pequeño (más cercano a cero que el menor número positivo representable).

Float32 boundaries:
  Maximum:  3.4028235e+38
  Minimum positive (normal): 1.175e-38
  Minimum positive (denorm): 1.401e-45
  Overflow:  anything > 3.4e38 becomes inf
  Underflow: anything < 1.4e-45 becomes 0.0

La función exp() es la principal fuente de overflow en ML:

exp(88.7)  = 3.40e+38   (barely fits in float32)
exp(89.0)  = inf         (overflow)
exp(-87.3) = 1.18e-38   (barely above underflow)
exp(-104)  = 0.0         (underflow to zero)

La función log() golpea en la dirección opuesta:

log(0.0)   = -inf
log(-1.0)  = nan
log(1e-45) = -103.3      (fine)
log(1e-46) = -inf        (input underflowed to 0, then log(0) = -inf)

En ML, exp() aparece en softmax, sigmoid y computaciones de probabilidad. log() aparece en entropía cruzada, log-verosimilitudes y divergencia KL. La combinación log(exp(x)) es un campo minado sin los trucos correctos.

El Truco del Log-Sum-Exp

Calcular log(sum(exp(x_i))) directamente es numéricamente peligroso. Si cualquier x_i es grande, exp(x_i) desborda. Si todos los x_i son muy negativos, cada exp(x_i) sufre underflow a cero y log(0) es -inf.

El truco: resta el valor máximo antes de exponenciar.

log(sum(exp(x_i))) = max(x) + log(sum(exp(x_i - max(x))))

Por que funciona: después de restar max(x), el mayor exponente es exp(0) = 1. Ningun overflow es posible. Al menos un término en la suma es 1, así que la suma es al menos 1, y log(1) = 0. Ningun underflow a -inf es posible.

Prueba:

log(sum(exp(x_i)))
= log(sum(exp(x_i - c + c)))                    (add and subtract c)
= log(sum(exp(x_i - c) * exp(c)))               (exp(a+b) = exp(a)*exp(b))
= log(exp(c) * sum(exp(x_i - c)))               (factor out exp(c))
= c + log(sum(exp(x_i - c)))                    (log(a*b) = log(a) + log(b))

Haz c = max(x) y el overflow queda eliminado.

Este truco aparece en todas partes en ML:

• Normalización del softmax
• Computación de la pérdida de entropía cruzada
• Suma de log-probabilidades en modelos de secuencia
• Mezcla de Gaussianas
• Inferencia variacional

Por que el Softmax Necesita el Truco de la Resta del Máximo

El softmax convierte logits en probabilidades:

softmax(x_i) = exp(x_i) / sum(exp(x_j))

Sin el truco, logits de [100, 101, 102] causan overflow:

exp(100) = 2.69e43
exp(101) = 7.31e43
exp(102) = 1.99e44
sum      = 2.99e44

These overflow float32 (max ~3.4e38)? No, 2.69e43 < 3.4e38? Actually:
exp(88.7) is already at the float32 limit.
exp(100) = inf in float32.

Con el truco, resta max(x) = 102:

exp(100 - 102) = exp(-2) = 0.135
exp(101 - 102) = exp(-1) = 0.368
exp(102 - 102) = exp(0)  = 1.000
sum = 1.503

softmax = [0.090, 0.245, 0.665]

Las probabilidades son idénticas. La computación es segura. Esto no es una optimización. Es un requisito para la corrección.

NaN e Inf: Detección y Prevencion

nan (Not a Number) e inf (infinito) se propagan viralmente por la computación. Un nan en una actualización de gradiente vuelve nan el peso, lo que vuelve nan toda salida subsiguiente. El entrenamiento está muerto en un solo paso.

Como aparece el inf:

• exp() de un número positivo grande
• División por cero: 1.0 / 0.0
• overflow de float32 en acumulaciones

Como aparece el nan:

• 0.0 / 0.0
• inf - inf
• inf * 0
• sqrt() de un número negativo
• log() de un número negativo
• Cualquier aritmética que involucre un nan existente

Detección:

import math

math.isnan(x)       # True if x is nan
math.isinf(x)       # True if x is +inf or -inf
math.isfinite(x)    # True if x is neither nan nor inf

Estrategias de prevencion:

1. Limita (clamp) las entradas de exp(): exp(clamp(x, -80, 80))
2. Agrega epsilon a los denominadores: x / (y + 1e-8)
3. Agrega epsilon dentro de log(): log(x + 1e-8)
4. Usa implementaciones estables (log-sum-exp, softmax estable)
5. Recorte de gradiente (gradient clipping) para evitar la explosión de pesos
6. Verifica nan/inf después de cada forward pass durante la depuración

Verificación Numérica de Gradiente

Los gradientes analíticos (de la retropropagación) pueden tener bugs. La verificación numérica de gradiente los verifica calculando gradientes con diferencias finitas.

La fórmula de la diferencia centrada:

df/dx ~= (f(x + h) - f(x - h)) / (2h)

Esto tiene precisión O(h^2), mucho mejor que la diferencia hacia adelante (f(x+h) - f(x)) / h, que es solo O(h).

Eligiendo h: demasiado grande y la aproximación queda errada. Demasiado pequeño y la cancelación catastrófica destruye la respuesta. h = 1e-5 a 1e-7 es lo típico.

La verificación: calcula la diferencia relativa entre los gradientes analíticos y numéricos.

relative_error = |grad_analytical - grad_numerical| / max(|grad_analytical|, |grad_numerical|, 1e-8)

Reglas practicas:

• relative_error < 1e-7: perfecto, el gradiente es correcto
• relative_error < 1e-5: aceptable, probablemente correcto
• relative_error > 1e-3: algo esta mal
• relative_error > 1: el gradiente esta completamente mal

Siempre verifica gradientes al implementar una nueva capa o función de pérdida. PyTorch provee torch.autograd.gradcheck() para esto.

Entrenamiento en Precisión Mixta

Las GPUs modernas tienen hardware especializado (Tensor Cores) que calcula multiplicaciones de matrices en float16 de 2 a 8 veces más rápido que en float32. El entrenamiento en precisión mixta explota esto:

1. Maintain float32 master copy of weights
2. Forward pass in float16 (fast)
3. Compute loss in float32 (prevents overflow)
4. Backward pass in float16 (fast)
5. Scale gradients to float32
6. Update float32 master weights

El problema con el entrenamiento puro en float16: los gradientes suelen ser muy pequeños (1e-8 o menores). float16 sufre underflow de cualquier valor por debajo de ~6e-8 a cero. Tu modelo deja de aprender porque todas las actualizaciones de gradiente son cero.

La corrección es el escalado de la pérdida (loss scaling):

1. Multiply loss by a large scale factor (e.g., 1024)
2. Backward pass computes gradients of (loss * 1024)
3. All gradients are 1024x larger (pushed above float16 underflow)
4. Divide gradients by 1024 before updating weights
5. Net effect: same update, but no underflow

El escalado dinámico de la pérdida ajusta el factor de escala automáticamente. Comienza con un valor grande (65536). Si los gradientes desbordan a inf, divídelo a la mitad. Si pasan N pasos sin overflow, duplícalo.

bfloat16 vs float16: Por que bfloat16 Gana para el Entrenamiento

float16:   [1 sign] [5 exponent]  [10 mantissa]
bfloat16:  [1 sign] [8 exponent]  [7 mantissa]

float16 tiene más precisión (10 bits de mantisa vs 7) pero rango limitado (max ~65.504). bfloat16 tiene menos precisión pero el mismo rango que float32 (max ~3.4e38).

Para entrenar redes neuronales:

• Las activaciones y los logits superan regularmente 65.504 durante picos de entrenamiento. float16 desborda; bfloat16 lo maneja.
• El escalado de la pérdida es obligatorio con float16 pero generalmente innecesario con bfloat16, porque su rango cubre el espectro de magnitud de los gradientes.
• bfloat16 es un truncamiento simple de float32: descarta los 16 bits inferiores de la mantisa. La conversión es trivial y sin pérdida en el exponente.

float16 se prefiere para inferencia, donde los valores están acotados y la precisión importa más. bfloat16 se prefiere para entrenamiento, donde el rango importa más. Por eso las TPUs y las GPUs NVIDIA modernas (A100, H100) tienen soporte nativo de bfloat16.

Recorte de Gradiente (Gradient Clipping)

Los gradientes explosivos pasan cuando los gradientes crecen exponencialmente a través de muchas capas (común en RNNs, redes profundas y transformers). Un solo gradiente grande puede corromper todos los pesos en un paso.

Dos tipos de recorte:

Recorte por valor: limita cada elemento del gradiente de forma independiente.

grad = clamp(grad, -max_val, max_val)

Simple, pero puede cambiar la dirección del vector gradiente.

Recorte por norma: escala el vector gradiente entero de modo que su norma no supere un umbral.

if ||grad|| > max_norm:
    grad = grad * (max_norm / ||grad||)

Preserva la dirección del gradiente. Es lo que hace torch.nn.utils.clip_grad_norm_(). Es la elección estándar.

Valores típicos: max_norm=1.0 para transformers, max_norm=0.5 para RL, max_norm=5.0 para redes más simples.

El recorte de gradiente no es un truco sucio. Es un mecanismo de seguridad. Sin el, un solo batch atípico puede producir un gradiente lo suficientemente grande para arruinar semanas de entrenamiento.

Capas de Normalización como Estabilizadores Numéricos

La normalización por batch (batch normalization), la normalización por capa (layer normalization) y la normalización RMS suelen presentarse como regularizadores que ayudan al entrenamiento a converger. También son estabilizadores numéricos.

Sin normalización, las activaciones pueden crecer o encoger exponencialmente a través de las capas:

Layer 1: values in [0, 1]
Layer 5: values in [0, 100]
Layer 10: values in [0, 10,000]
Layer 50: values in [0, inf]

La normalización recentra y reescala las activaciones en cada capa:

LayerNorm(x) = (x - mean(x)) / (std(x) + epsilon) * gamma + beta

El epsilon (típicamente 1e-5) evita la división por cero cuando todas las activaciones son idénticas. Los parámetros aprendidos gamma y beta dejan que la red restaure cualquier escala que necesite.

Esto mantiene los valores en un rango numéricamente seguro a lo largo de la red, evitando tanto el overflow en el forward pass como la explosión de gradiente en el backward pass.

Bugs Numéricos Comunes de ML

Bug: La pérdida se vuelve NaN después de algunas épocas. Causa: los logits crecieron demasiado, el softmax desbordo. O la tasa de aprendizaje es demasiado alta y los pesos divergieron. Corrección: usa softmax estable (resta del máximo), reduce la tasa de aprendizaje, agrega recorte de gradiente.

Bug: La pérdida queda atascada en log(num_classes). Causa: las salidas del modelo son probabilidades casi uniformes. A menudo significa que los gradientes se están desvaneciendo o que el modelo no está aprendiendo nada. Corrección: verifica que las etiquetas de los datos sean correctas, verifica la función de pérdida, verifica ReLUs muertas.

Bug: La exactitud de validación esta 1-3% por debajo de lo esperado. Causa: precisión mixta sin el escalado adecuado de la pérdida. El underflow de gradiente pone en cero silenciosamente las actualizaciones pequeñas. Corrección: activa el escalado dinámico de la pérdida, o cambia a bfloat16.

Bug: Las normas de gradiente son 0.0 para algunas capas. Causa: neuronas ReLU muertas (todas las entradas negativas), o underflow de float16. Corrección: usa LeakyReLU o GELU, usa escalado de gradiente, verifica la inicialización de pesos.

Bug: El modelo funciona en una GPU pero da resultados diferentes en otra. Causa: orden de acumulación de punto flotante no determinista. Las reducciones paralelas de la GPU suman en órdenes diferentes en hardware diferente, y la suma de punto flotante no es asociativa. Corrección: acepta pequeñas diferencias (1e-6), o define torch.use_deterministic_algorithms(True) y acepta la penalización de velocidad.

Bug: exp() devuelve inf en la computación de la pérdida. Causa: logits crudos pasados a exp() sin el truco de la resta del máximo. Corrección: usa torch.nn.functional.log_softmax(), que implementa log-sum-exp internamente.

Bug: El entrenamiento diverge después de cambiar de float32 a float16. Causa: float16 no puede representar magnitudes de gradiente por debajo de 6e-8 ni activaciones por encima de 65.504. Corrección: usa precisión mixta con escalado de la pérdida (AMP), o usa bfloat16.

Construye

Paso 1: Demuestra los límites de precisión del punto flotante

print("=== Floating Point Precision ===")
print(f"0.1 + 0.2 = {0.1 + 0.2}")
print(f"0.1 + 0.2 == 0.3? {0.1 + 0.2 == 0.3}")
print(f"Difference: {(0.1 + 0.2) - 0.3:.2e}")

Paso 2: Implementa softmax ingenuo vs estable

import math

def softmax_naive(logits):
    exps = [math.exp(z) for z in logits]
    total = sum(exps)
    return [e / total for e in exps]

def softmax_stable(logits):
    max_logit = max(logits)
    exps = [math.exp(z - max_logit) for z in logits]
    total = sum(exps)
    return [e / total for e in exps]

safe_logits = [2.0, 1.0, 0.1]
print(f"Naive:  {softmax_naive(safe_logits)}")
print(f"Stable: {softmax_stable(safe_logits)}")

dangerous_logits = [100.0, 101.0, 102.0]
print(f"Stable: {softmax_stable(dangerous_logits)}")
# softmax_naive(dangerous_logits) would return [nan, nan, nan]

Paso 3: Implementa log-sum-exp estable

def logsumexp_naive(values):
    return math.log(sum(math.exp(v) for v in values))

def logsumexp_stable(values):
    c = max(values)
    return c + math.log(sum(math.exp(v - c) for v in values))

safe = [1.0, 2.0, 3.0]
print(f"Naive:  {logsumexp_naive(safe):.6f}")
print(f"Stable: {logsumexp_stable(safe):.6f}")

large = [500.0, 501.0, 502.0]
print(f"Stable: {logsumexp_stable(large):.6f}")
# logsumexp_naive(large) returns inf

Paso 4: Implementa entropía cruzada estable

def cross_entropy_naive(true_class, logits):
    probs = softmax_naive(logits)
    return -math.log(probs[true_class])

def cross_entropy_stable(true_class, logits):
    max_logit = max(logits)
    shifted = [z - max_logit for z in logits]
    log_sum_exp = math.log(sum(math.exp(s) for s in shifted))
    log_prob = shifted[true_class] - log_sum_exp
    return -log_prob

logits = [2.0, 5.0, 1.0]
true_class = 1
print(f"Naive:  {cross_entropy_naive(true_class, logits):.6f}")
print(f"Stable: {cross_entropy_stable(true_class, logits):.6f}")

Paso 5: Verificación de gradiente

def numerical_gradient(f, x, h=1e-5):
    grad = []
    for i in range(len(x)):
        x_plus = x[:]
        x_minus = x[:]
        x_plus[i] += h
        x_minus[i] -= h
        grad.append((f(x_plus) - f(x_minus)) / (2 * h))
    return grad

def check_gradient(analytical, numerical, tolerance=1e-5):
    for i, (a, n) in enumerate(zip(analytical, numerical)):
        denom = max(abs(a), abs(n), 1e-8)
        rel_error = abs(a - n) / denom
        status = "OK" if rel_error < tolerance else "FAIL"
        print(f"  param {i}: analytical={a:.8f} numerical={n:.8f} "
              f"rel_error={rel_error:.2e} [{status}]")

def f(params):
    x, y = params
    return x**2 + 3*x*y + y**3

def f_grad(params):
    x, y = params
    return [2*x + 3*y, 3*x + 3*y**2]

point = [2.0, 1.0]
analytical = f_grad(point)
numerical = numerical_gradient(f, point)
check_gradient(analytical, numerical)

Usalo

Simulación de precisión mixta

import struct

def float32_to_float16_round(x):
    packed = struct.pack('f', x)
    f32 = struct.unpack('f', packed)[0]
    packed16 = struct.pack('e', f32)
    return struct.unpack('e', packed16)[0]

def simulate_bfloat16(x):
    packed = struct.pack('f', x)
    as_int = int.from_bytes(packed, 'little')
    truncated = as_int & 0xFFFF0000
    repacked = truncated.to_bytes(4, 'little')
    return struct.unpack('f', repacked)[0]

Recorte de gradiente

def clip_by_norm(gradients, max_norm):
    total_norm = math.sqrt(sum(g**2 for g in gradients))
    if total_norm > max_norm:
        scale = max_norm / total_norm
        return [g * scale for g in gradients]
    return gradients

grads = [10.0, 20.0, 30.0]
clipped = clip_by_norm(grads, max_norm=5.0)
print(f"Original norm: {math.sqrt(sum(g**2 for g in grads)):.2f}")
print(f"Clipped norm:  {math.sqrt(sum(g**2 for g in clipped)):.2f}")
print(f"Direction preserved: {[c/clipped[0] for c in clipped]} == {[g/grads[0] for g in grads]}")

Detección de NaN/Inf

def check_tensor(name, values):
    has_nan = any(math.isnan(v) for v in values)
    has_inf = any(math.isinf(v) for v in values)
    if has_nan or has_inf:
        print(f"WARNING {name}: nan={has_nan} inf={has_inf}")
        return False
    return True

check_tensor("good", [1.0, 2.0, 3.0])
check_tensor("bad",  [1.0, float('nan'), 3.0])
check_tensor("ugly", [1.0, float('inf'), 3.0])

Consulta code/numerical.py para implementaciones completas con todos los casos límite demostrados.

Entregalo

Esta lección produce:

• code/numerical.py con softmax estable, log-sum-exp, entropía cruzada, verificación de gradiente y simulación de precisión mixta
• outputs/prompt-numerical-debugger.md para diagnosticar NaN/Inf y problemas numéricos en el entrenamiento

Estas implementaciones estables reaparecen en la Fase 3 al construir el bucle de entrenamiento y en la Fase 4 al implementar mecanismos de atención.

Ejercicios

1. Cancelación catastrófica. Calcula la varianza de [1000000.0, 1000001.0, 1000002.0] usando la fórmula ingenua E[x^2] - E[x]^2 en float32. Luego calcúlala usando el algoritmo online de Welford. Compara los errores contra la varianza verdadera (0.6667).

2. Caza de precisión. Encuentra el menor valor positivo de float32 x tal que 1.0 + x == 1.0 en Python. Ese es el epsilon de máquina. Verifica que coincida con numpy.finfo(numpy.float32).eps.

3. Casos límite del log-sum-exp. Prueba tu función logsumexp_stable con: (a) todos los valores iguales, (b) un valor mucho mayor que los demas, (c) todos los valores muy negativos (-1000). Verifica que de resultados correctos donde la versión ingenua falla.

4. Verificación de gradiente en una capa de red neuronal. Implementa una sola capa lineal y = Wx + b y su backward pass analítico. Usa numerical_gradient para verificar la corrección para una matriz de pesos 3x2.

5. Experimento de escalado de la pérdida. Simula entrenamiento con float16: crea gradientes aleatorios en el rango [1e-9, 1e-3], conviértelos a float16 y mide que fracción se vuelve cero. Luego aplica el escalado de la pérdida (multiplica por 1024), convierte a float16, reescala de vuelta y mide la fracción de ceros de nuevo.

Términos Clave

Lectura Adicional

• What Every Computer Scientist Should Know About Floating-Point Arithmetic (Goldberg 1991) -- la referencia definitiva, densa pero completa
• Mixed Precisión Training (Micikevicius et al., 2018) -- el paper de NVIDIA que introdujo el escalado de la pérdida para el entrenamiento en float16
• AMP: Automatic Mixed Precisión (PyTorch docs) -- guia practica de precisión mixta en PyTorch
• bfloat16 format (Google Cloud TPU docs) -- por que Google eligió este formato para las TPUs
• Kahan Summation (Wikipedia) -- algoritmo para reducir el error de redondeo en sumas de punto flotante