Skip to content

English · Español

Lab 02 — Medir el error de suma y arreglarlo con Kahan

Objetivo: ver una suma fp32 ingenua dejar de crecer, luego ver a la suma de Kahan salvarla. Ancla: sumar las 600 probabilidades por forma del vocabulario verbal de §A13 y verificar que el total iguala 1.0 con precisión ε.

Tiempo estimado: 60–90 minutos.

Prerrequisito: teoría 03-summation-and-cancellation.md leída.

Lo que produces

Un directorio experiments/02-kahan-vs-naive/ que contiene:

  • summation.py — implementaciones de naive_sum, pairwise_sum (emulación del predeterminado de NumPy), kahan_sum, neumaier_sum.
  • experiment.py — driver que ejecuta los cuatro sobre las tres cargas de trabajo de test (abajo) y registra el error frente a una referencia fp64.
  • results.json — números de error por carga de trabajo.
  • error_growth.png — gráfico log-log del error acumulado frente a N.
  • manifest.json.
  • README.md — interpretación.

Las tres cargas de trabajo

Carga 1 — probabilidades uniformes del vocabulario de §A13

Un array fp32 de longitud 600, todas las entradas iguales a 1.0 / 600.0. La suma debería ser exactamente 1.0.

probs = np.full(600, 1.0/600.0, dtype=np.float32)

Suma ingenua, suma por pares, suma de Kahan. Registra cada resultado y el error relativo a 1.0.

Carga 2 — 10⁶ muestras de una distribución de magnitud pequeña

Genera N = 10⁶ valores fp32, cada uno ~1/600 + δ_i donde δ_i ~ Normal(0, 1e-5). Usa np.random.default_rng(42).

N = 1_000_000
rng = np.random.default_rng(42)
x = (1.0/600.0 + rng.standard_normal(N).astype(np.float32) * 1e-5).astype(np.float32)
expected = x.astype(np.float64).sum()  # fp64 reference

Ejecuta las cuatro estrategias de suma. Grafica el error vs N en incrementos (barre N de 10³ a 10⁶ por factores de 10) — esto es error_growth.png.

Carga 3 — uno grande + muchos diminutos

Construye una suma de magnitud mixta en el peor caso: un 1.0e6 seguido por 10⁶ unos.

x = np.array([1.0e6] + [1.0] * 1_000_000, dtype=np.float32)
expected = 1_001_000.0  # exact in fp32 (representable)

Predice: la suma ingenua pierde cada 1.0 después de que el total corriente exceda 2²³ ≈ 8.4e6 (el punto donde 1.0 es más pequeño que el ULP del total corriente). Kahan los mantiene todos. Verifica.

TODOs

Bloque A — implementaciones

summation.py:

def naive_sum(arr):
    total = 0.0
    for x in arr:
        total += x
    return total

def pairwise_sum(arr):
    # Recursive halving — like np.sum's default
    if len(arr) <= 8:
        return float(arr.sum())   # base case OK to use np.sum here
    mid = len(arr) // 2
    return pairwise_sum(arr[:mid]) + pairwise_sum(arr[mid:])

def kahan_sum(arr):
    # TODO: implement per theory/03 §"the fix"
    ...

def neumaier_sum(arr):
    # TODO: implement per theory/03 §"Neumaier's variant"
    ...

Todas deben aceptar un array 1D np.float32 y devolver un float (el fp64 nativo de Python está bien para el tipo de retorno — la disciplina es interna).

Variante: escribe también fp64_accumulator_naive(arr) que haga total = np.float64(0.0) y acumule en fp64. Esta es la comparación de "la solución correcta más barata".

Bloque B — predecir antes de medir

En README.md, antes de ejecutar, predice para cada carga de trabajo:

  • Magnitud del error ingenuo fp32.
  • Magnitud del error por pares fp32.
  • Magnitud del error Kahan fp32.
  • Magnitud del error de acumulador-fp64-sobre-entrada-fp32.

Usa la teoría 03 § "La solución" para predecir. Luego ejecuta y compara.

Bloque C — carga 1: vocabulario de §A13

Para el vector de probabilidades uniformes de longitud 600, los cuatro métodos deberían dar errores cómodamente por debajo de ε. Declara cuál método es mejor y por qué.

🇪🇸 La intuición: 600 sumas de magnitudes idénticas tienen error naïve O(600 · ε) ≈ 7e-5. Kahan: O(ε) ≈ 1e-7. fp64 accumulator: O(ε_{fp64}) ≈ 2e-16. Si crees que total debe ser bit-exacto 1.0, ninguno de los métodos te lo garantiza — los δ_i de redondeo en 1/600 mismo ya rompen la igualdad exacta.

Bloque D — carga 2: error vs N

Genera el array de 10⁶. Para N ∈ [10³, 10⁴, 10⁵, 10⁶], calcula el error de cada método de suma y regístralo. Grafica error_growth.png:

  • eje x: N (log).
  • eje y: error absoluto frente a la referencia fp64 (log).
  • Una curva por método.
  • Anotaciones: cotas teóricas O(N · ε) para ingenuo, O(log N · ε) para por pares, O(ε) para Kahan.

La forma esperada es ingenuo creciendo linealmente, por pares como log N, Kahan plano. Comenta si tus medidas coinciden en README.md.

Bloque E — carga 3: caso catastrófico

Ejecuta los cuatro métodos sobre el array [1e6] + [1.0] * 10⁶. Predice y verifica:

  • naive_sum devuelve... ¿qué? ¿En qué índice se vuelve un no-op sumar 1.0 al total corriente?
  • pairwise_sum es algo mejor pero todavía no perfecto.
  • kahan_sum y neumaier_sum deberían recuperar ambos la respuesta exacta (o dentro de 1 ULP).

Calcula el índice en el que ingenuo deja de crecer: cuando total > 2²³ ≈ 8.4e6. Para nuestra carga, total alcanza 1e6 + i × 1.0, así que ingenuo deja de crecer en i ≈ 7.4e6 - 1e6 = ... — pero el array solo tiene 1e6 de longitud, así que para este caso específico el asunto es más sutil. Resuélvelo, escribe la respuesta en README.md.

Bloque F — la recomendación práctica

En README.md, concluye con un párrafo: cuál de estos métodos usaría Borja en código de entrenamiento real, cuándo. Anclate en este escenario:

Entrenando MiniGPT sobre el vocabulario de §A13, acumulas pérdidas por batch sobre 100,000 batches × 64 ejemplos cada uno = 6.4 millones de sumandos por época. Cada sumando es una pérdida fp32 positiva en [0, 10]. La acumulación ingenua da ~5% de error relativo. ¿Qué solución?

Respuesta esperada: total = np.float64(0.0) y acumular en fp64. Cambio de código trivial, sin complejidad de Kahan, precisión suficiente. Recurre a Kahan solo si estás atascado en fp32 (algunos caminos de inferencia embebida).

Restricciones

  • Nada de np.sum en naive_sum, kahan_sum o neumaier_sum. Hazte el bucle. El punto es ver la acumulación.
  • Los bucles puros de Python son lentos. Para N = 10⁶ pueden tardar minutos. Para el gráfico de la carga 2 puedes tolerar la lentitud o implementar los bucles con numba.njit (opcional; fíjalo solo en el grupo experiments si vas por este camino — no hinches el entorno).
  • Usa siempre la referencia fp64. Nunca la "respuesta matemática esperada" — la suma fp32 tiene un suelo genuino de error de redondeo que fp64 revela.

Condiciones de parada

Hecho cuando:

  1. Las cuatro implementaciones de suma funcionan y se auto-testean (la suma de [1.0]*10 devuelve 10.0).
  2. results.json contiene errores para las tres cargas de trabajo × cuatro métodos.
  3. error_growth.png muestra el escalado predicho.
  4. README.md hace la recomendación práctica explícitamente.
  5. Puedes decir, en una frase, cuándo Kahan vale su complejidad frente a un acumulador fp64.

Escollos

  • Mezclar fp32 y fp64. 1.0/600.0 en Python es fp64. np.float32(1.0/600.0) trunca a fp32. Los arrays deben ser np.float32 para que el experimento demuestre el error fp32.
  • Llamar arr.sum() accidentalmente. Los arrays np.float32 usan por defecto la suma por pares de NumPy — mucho mejor que ingenuo. El punto es el bucle ingenuo. No tomes atajos.
  • Varianza de tiempo. Si implementas naive_sum en Python puro, será 100-1000× más lento que np.sum. Eso es esperable; el rendimiento no es la medición aquí.
  • Suelo de referencia fp64. Incluso la referencia fp64 tiene ~10⁻¹⁶ de error relativo. Para la carga 2 tu error fp32 Kahan de ~10⁻⁷ está bien, pero no esperes que coincida con fp64 por identidad de bits.

Cuándo consultar solutions/

Tras commitear los cinco archivos. Solución en solutions/02-summation-ref.md (escrita al abrir la fase).

Siguiente laboratorio: lab/03-quantization-preview.md.