NumPy como escape hatch para cálculos numéricos

Cuando tienes un cuello de botella en cálculos numéricos, la primera pregunta no es “¿cómo optimizo este loop?”, sino “¿por qué estoy usando un loop de Python aquí?”. NumPy existe precisamente para que esa pregunta tenga respuesta.

Qué es realmente NumPy (y por qué es rápido)

Un array de NumPy (ndarray) no es una lista de Python con superpoderes. Es un bloque contiguo de memoria homogénea —todos los elementos del mismo tipo C— con metadatos de forma y stride encima. Cuando llamas a arr * 2, Python no itera: delega esa multiplicación a una función compilada en C que opera directamente sobre ese bloque de memoria en un solo dispatch. El bytecode de Python solo ejecuta una instrucción; el bucle real ocurre en C, sin overhead de referencia a objetos, sin boxing/unboxing, sin GIL por elemento.

Eso explica el salto de rendimiento. Con una lista de Python, [x * 2 for x in data] ejecuta un ciclo completo del intérprete por cada elemento: resolución de nombre, creación de objeto int, gestión de refcounts. Con un ndarray, todo eso desaparece.

Broadcasting es el mecanismo que permite operar sobre arrays de formas distintas sin copiar datos ni escribir loops. NumPy extiende dimensiones virtuales siguiendo reglas específicas: dos arrays son compatibles eje por eje si sus tamaños son iguales o uno de ellos es 1. Esto no es azúcar sintáctica: es la forma correcta de expresar esas operaciones, y además es la más rápida.

Un loop explícito sobre un ndarrayfor x in arr— es siempre una anti-señal. Indica que estás usando NumPy como almacenamiento pero ejecutando la operación en Python. Obtienes lo peor de ambos mundos: la rigidez del array sin el beneficio de la velocidad C. Cuando te encuentres escribiendo ese loop, la pregunta es: ¿qué operación de array expresa esto directamente?

Ejemplo completo

import numpy as np
import time

# Simulación: normalización de un lote de imágenes en escala de grises.
# Shape: (N, H, W) → N imágenes de H×W píxeles, valores uint8 [0, 255].
rng = np.random.default_rng(42)
images = rng.integers(0, 256, size=(1000, 128, 128), dtype=np.uint8)

# ── Versión Python puro ───────────────────────────────────────────────
def normalize_python(imgs):
    result = []
    for img in imgs:
        row_result = []
        for row in img:
            row_result.append([px / 255.0 for px in row])
        result.append(row_result)
    return result

# ── Versión NumPy ─────────────────────────────────────────────────────
def normalize_numpy(imgs):
    # El cast a float32 antes de dividir evita que NumPy promueva a float64,
    # que consumiría el doble de memoria sin ganancia real para redes neuronales.
    return imgs.astype(np.float32) / 255.0

# ── Broadcasting: restar la media por canal ───────────────────────────
# images_f: (1000, 128, 128). mean_per_image: (1000, 1, 1).
# NumPy estira las dimensiones de size=1 para hacer la resta elemento a elemento
# sin materializar ninguna copia intermedia del tamaño completo.
images_f = images.astype(np.float32) / 255.0
mean_per_image = images_f.mean(axis=(1, 2), keepdims=True)  # shape: (1000, 1, 1)
centered = images_f - mean_per_image                         # broadcasting aquí

# ── Benchmark ────────────────────────────────────────────────────────
t0 = time.perf_counter()
normalize_python(images)
t_python = time.perf_counter() - t0

t0 = time.perf_counter()
normalize_numpy(images)
t_numpy = time.perf_counter() - t0

print(f"Python puro: {t_python:.3f}s")
print(f"NumPy:       {t_numpy:.4f}s")
print(f"Speedup:     {t_python / t_numpy:.1f}×")

# ── Verificación de coherencia del broadcasting ───────────────────────
print(f"\nForma original:        {images_f.shape}")
print(f"Media por imagen:      {mean_per_image.shape}")
print(f"Resultado centrado:    {centered.shape}")
assert centered.shape == images_f.shape

Qué está pasando en cada decisión

astype(np.float32) no es un detalle menor. uint8 no puede representar fracciones, así que sin ese cast NumPy truncaría los resultados. Usar float32 en lugar de float64 —el default de NumPy cuando ve un literal 1.0— reduce el consumo de memoria a la mitad. Para cálculos de ML, float32 es el estándar; para análisis estadístico de alta precisión, quizás necesitas float64. La decisión es tuya, pero tienes que tomarla conscientemente.

El eje del broadcasting merece atención. mean(axis=(1, 2), keepdims=True) colapsa las dimensiones de altura y anchura pero mantiene la estructura de 3 ejes con shape (1000, 1, 1). Sin keepdims=True obtendrías shape (1000,), y la resta images_f - mean_per_image fallaría o haría algo inesperado porque NumPy alinea dimensiones desde la derecha. keepdims es la forma explícita de comunicar tu intención al motor de broadcasting.

La división / 255.0 aplica la operación a todos los millones de elementos en C. No hay loop Python visible, no hay objeto float de Python creado por elemento. El resultado es un nuevo array contiguo en memoria, listo para ser pasado a SciPy, pandas o PyTorch, que entienden ndarray de forma nativa porque construyen sobre la misma infraestructura de buffer.

Errores que debes conocer

Error: Iterar sobre un array para aplicar una transformación elemento a elemento, perdiendo toda ventaja de NumPy.

# ❌ Wrong
result = np.array([np.sqrt(x) for x in arr])

# ✅ Right
result = np.sqrt(arr)

np.sqrt —como todas las ufuncs de NumPy— opera sobre el array completo en C. La list comprehension crea un objeto Python por elemento antes de construir el array de vuelta.

Error: Asumir que el broadcasting ocurrirá cuando las formas son incompatibles, obteniendo un error críptico o, peor, un resultado silenciosamente incorrecto.

# ❌ Wrong — shape (1000, 128, 128) - shape (1000,) → error o broadcast incorrecto
mean_wrong = images_f.mean(axis=(1, 2))          # shape: (1000,)
centered_wrong = images_f - mean_wrong            # ValueError: no puede broadcast

# ✅ Right — keepdims preserva los ejes para que el broadcasting sea inequívoco
mean_right = images_f.mean(axis=(1, 2), keepdims=True)  # shape: (1000, 1, 1)
centered_right = images_f - mean_right

NumPy alinea dimensiones desde el eje más a la derecha. (1000, 128, 128) contra (1000,) intenta emparejar 128 con 1000, que no son compatibles. keepdims=True produce (1000, 1, 1), que sí encaja sin ambigüedad.

Error: Encadenar operaciones de array sin considerar que cada una materializa un array intermedio completo en memoria.

# ❌ Wrong — crea tres arrays temporales del tamaño completo
result = (arr - arr.mean()) / arr.std() * scale + offset

# ✅ Right — cuando el tamaño importa, usa out= o numexpr para fusionar operaciones
# Para arrays grandes, considera numexpr o escribir la expresión en una sola ufunc.
result = np.empty_like(arr, dtype=np.float32)
np.subtract(arr, arr.mean(), out=result)
result /= arr.std()

Con arrays pequeños esto es irrelevante. Con arrays de varios gigabytes, cada operación intermedia puede provocar presión de memoria o incluso un OOM. El parámetro out= de las ufuncs permite operar in-place y evitar esas copias.

168

Dejar un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Scroll al inicio