Estabilidad numérica

01 · Punto de partida

Misma matemática, distinto programa

Los Temas 39–42 te dieron el diagnóstico: grilla finita (39), ½ ulp por operación (40), sesgo al cortar (41) y suma que acumula (42). Ahora la cura: dos programas que en el pizarrón son idénticos pueden comportarse opuesto en la máquina. El ejemplo estrella: la varianza “de libro” Σx²/N − media² resta dos gigantes casi iguales y devuelve varianzas negativas; el incremental de Welford actualiza media y dispersión paso a paso y queda clavado en la verdad.

La lección general: el error de redondeo es inevitable, pero amplificarlo es opcional. Un algoritmo estable lo deja en 1–2 ulps; uno inestable lo multiplica por el número de condición hasta volverlo basura.

• ¿Qué distingue un problema mal condicionado de un algoritmo inestable?
• ¿Por qué restar casi-iguales (cancelación) es la放大器 del error?
• ¿Cómo se calcula varianza sin restar gigantes?
• ¿Qué recetas estables uso en simulación: Welford, dos pasadas, log?

02 · Definición

Condición del problema, estabilidad del algoritmo

03 · Anatomía

Las tres recetas estables

Casi toda inestabilidad es una resta de casi-iguales disfrazada. Tres formas de desarmarla:

Centrar antes de operar (dos pasadas)

Primero la media, después Σ(x − media)²/N: los números chicos se restan entre sí, no dos gigantes. Dos recorridas, error mínimo.

Acumular incremental (Welford)

Media y M₂ se actualizan dato a dato sin guardar el vector: una pasada, estable y online. El estándar para streams de simulación.

Reformular la expresión

Cambiar la fórmula por una equivalente sin cancelación: 1 − cos(x) → 2·sin²(x/2), log(1+x) con log1p, productos con fma.

04 · Ejemplos

Tres inestabilidades clásicas

1. 1
   Varianza negativa.

   Datos 1e8 ± 2: la de libro resta 1,0000000000000004e16 − 1,0000000000000000e16 y puede dar −4,0. Welford da 2,0 clavado.

2. 2
   1 − cos(1e−8).

   Ingenuo: 1 − 0,9999999999999999 = 0,0 (todo cancelado). Estable: 2·sin²(5e−9) ≈ 5e−17 correcto.

3. 3
   Promedio de 1 M de energías.

   Sumar todo y dividir arrastra absorción (Tema 42); el promedio incremental (μ += (x−μ)/n) nunca forma el gigante y es estable por construcción.

05 · Demostración en Python

La varianza que da negativa

El mismo vector, dos funciones: la de libro y Welford. Con media 1e8 la primera se rompe y la segunda ni se inmuta.

def var_libro(xs):
    n = len(xs)
    return sum(x * x for x in xs) / n - (sum(xs) / n) ** 2

def var_welford(xs):
    n = mu = m2 = 0
    for x in xs:
        n += 1
        d = x - mu
        mu += d / n
        m2 += d * (x - mu)
    return m2 / n


base, vals = 1e8, [1e8 + d for d in [-2, -1, 0, 1, 2]] * 400
print("libro:  ", var_libro(vals))
print("welford:", var_welford(vals))

Libro ≈ negativa o disparatada; Welford ≈ 2,0: la resta de gigantes delata al inestable.

Promedio incremental (bonus estable)

mu = 0.0
for n, x in enumerate(vals, start=1):
    mu += (x - mu) / n
print(mu)

06 · Exploración

Laboratorio: la varianza que explota con la media

Datos con varianza verdadera 2,0 (patrón −2…+2 repetido) montados sobre una media base M. La curva naranja es el error de la fórmula de libro y la verde el de Welford, en función de M (eje log). Elegí M y N y mirá dónde se rompe la ingenua.

El eje x es la media base (log), el eje y el error absoluto en log. La ingenua sube en diagonal con M²; la estable queda abajo pegada al cero.

Preguntas para explorar

1. Con M = 0: ¿ambas dan 2,0? ¿Por qué la cancelación no muerde ahí?
2. Con M = 1e8: ¿la de libro da negativa, cero o gigante? ¿Qué dígitos se cancelaron?
3. Subí N de 500 a 5000 con M = 1e8: ¿la ingenua mejora o empeora? ¿Y Welford?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: romper la de libro

Reproducí la varianza negativa con base 1e8 y confirmá que Welford da 2,0.

vals = [1e8 + d for d in [-2, -1, 0, 1, 2]] * 400
n = len(vals)
print(sum(x * x for x in vals) / n - (sum(vals) / n) ** 2)

Ejercicio 2: 1 − cos(x) que da cero

Mostrá que la forma ingenua colapsa con x chico y que 2·sin²(x/2) la salva.

import math

x = 1e-8
print(1 - math.cos(x))
print(2 * math.sin(x / 2) ** 2)

Ejercicio 3: Welford desde cero

Escribí welford que devuelva media y varianza en una pasada y probalo con base 1e6.

def welford(xs):
    n = mu = m2 = 0
    for x in xs:
        n += 1
        d = x - mu
        mu += d / n
        m2 += d * (x - mu)
    return mu, m2 / n


vals = [1e6 + d for d in [-2, -1, 0, 1, 2]] * 200
print(welford(vals))

09 · Síntesis

Ideas para recordar

En el próximo tema volvemos al azar: cómo la calidad del generador sesga todos los resultados de la simulación.