Skip to content

English · Español

Lab 01 — Rendimiento de matmul: ingenuo vs por bloques vs einsum vs BLAS

Objetivo: medir el gap 50× entre el matmul triple-bucle de Python y np.matmul, y explicarlo como un hecho de roofline (Fase 1) compuesto por vectorización + evitar el intérprete.

Tiempo estimado: 90–120 minutos.

Prerrequisito: theory/02-matmul-and-shapes.md, lab/00-shapes-by-hand.md completos, resultado del roofline de la Fase 1 accesible desde experiments/01-*/.

Lo que produces

Un directorio experiments/03-matmul-perf/ que contiene:

  • bench.py — tu script de benchmark.
  • results.json — mediciones a varios N.
  • perf.png — gráfico log-log de GFLOPS vs N, cuatro curvas + overlay del roofline.
  • manifest.json.
  • README.md — interpretación, 2-3 párrafos.

No se introduce ningún módulo src/ en esta fase. Las cuatro variantes de matmul viven como funciones planas en bench.py (o un archivo helper hermano dentro del directorio del experimento). La implementación se promueve a src/minigrad/linalg.py en la Fase 7 cuando llega el consumidor de autograd — hasta entonces, el estado de script scratch es lo correcto.

Parte A — Implementar las cuatro variantes de matmul

Dentro de experiments/03-matmul-perf/:

  • matmul_naive(A, B) — triple-bucle Python puro (con arrays NumPy como contenedores, sin listas Python; pero el multiply-add escalar interior se interpreta). FP32 en todo.
  • matmul_blocked(A, B, block_size=64) — versión por bloques de seis bucles. Documenta el tamaño de bloque como parámetro.
  • matmul_einsum(A, B) — una sola llamada a np.einsum('ik,kj->ij', A, B).
  • np.matmul(A, B) — la referencia.
  • Verifica que los cuatro coinciden con A @ B con rtol=1e-5, atol=1e-6 sobre entradas aleatorias (M=K=N=64). (Las diferencias de acumulador en matmul son reales pero acotadas.)

Restricciones

  • Pure NumPy. Sin scipy.linalg, sin torch, sin numba, sin cython.
  • Sin vectorización temprana en matmul_naive. Si escribes C[i, k] += A[i, k] * B[k, j] dentro de tres bucles for, lo has hecho bien. Si escribes C[i] += A[i, k] * B[k] (haciendo broadcast por j), accidentalmente has vectorizado el bucle interno — los GFLOPS medidos serán demasiado altos y el lab pierde su sentido.
  • Determinismo. np.random.default_rng(42).standard_normal(...).astype(np.float32) para las entradas.
  • bench.py también escribe manifest.json con seed + versiones + governor de CPU + ajustes de threads.

Parte B — El benchmark

bench.py mide el tiempo de pared de cada método a varios N:

  • Tamaños: N ∈ {64, 128, 256, 512, 1024}. Matrices cuadradas. fp32.
  • Métodos: naive, blocked, einsum, np.matmul. Opcionalmente también blocked_32, blocked_128 para curiosidad de tuning.
  • Iteraciones: suficientes para tragarse el ruido del timer — mínimo 5 reps, con el tiempo total transcurrido por (método, N) de al menos 1 segundo. Para naive en N=1024 puede que necesites ejecutar menos iteraciones (o saltar el tamaño mayor y extrapolar desde N=512 por T_naive(N) ≈ T_naive(N/2) × 8 ya que el coste es O(N³)).
  • Warm-up: una iteración no cronometrada por (método, N) antes de medir.
  • GFLOPS: (2 · N³) / elapsed_seconds / 1e9. Guarda como results.json.

Cota de cordura

Deberías ver aproximadamente:

Método N=128 N=1024
naive ~3 MFLOPS ~3 MFLOPS (el intérprete Python domina independientemente de N)
blocked ~3-5 MFLOPS similar (sigue siendo el intérprete Python)
einsum ~10-30 GFLOPS ~20-60 GFLOPS
np.matmul ~20-50 GFLOPS ~50-100 GFLOPS

Si naive corre a >100 MFLOPS, vectorizaste accidentalmente — vuelve a revisar que no haya ops broadcast de NumPy en el bucle interior. Si np.matmul corre a <10 GFLOPS, OpenBLAS no está enlazado o estás throttled. Comprueba np.show_config() y el governor de CPU.

Parte C — Gráfico

perf.png:

  • eje x: N (escala log, ticks enteros en 64, 128, 256, 512, 1024).
  • eje y: GFLOPS (escala log).
  • Cuatro curvas, una por método.
  • Overlay: línea horizontal en los GFLOPS pico fp32 medidos del i5-8250U del experimento de roofline de la Fase 1. Calcula el cruce con el techo de memoria. Anota.
  • Título: "Rendimiento de matmul en i5-8250U, fp32, single-threaded."

Parte D — Interpretación del README

En README.md, responde:

  1. ¿Cuál es el speedup de np.matmul sobre matmul_naive en N=1024? Indícalo como multiplicador.
  2. Descompón el speedup en "magia BLAS" vs "sin intérprete Python". Pista: einsum típicamente captura la mayor parte del speedup de BLAS (despacha a las mismas rutinas), pero np.matmul puede adelantarse en N pequeño por diferencias en pases del optimizador. Así que aproximadamente: eliminar-intérprete ≈ 100×, SIMD ≈ 8×, cache blocking ≈ 5-10×, multi-thread (apagado aquí) sería otro 4×.
  3. ¿Dónde se sitúa np.matmul en el roofline de la Fase 1? Debería caer cerca o encima del techo de memoria en N=1024 (porque BLAS reusa la cache agresivamente, elevando la intensidad aritmética efectiva bien por encima del piso ingenuo 0.25 FLOP/byte). Indica la intensidad aritmética que alcanza.
  4. ¿matmul_blocked supera a matmul_naive en Python? Si no, ¿por qué? (Pista: la sobrecarga del intérprete Python eclipsa el coste de cache miss que el bloqueo pretende amortizar. El bloqueo es una optimización C/asm; en Python puro es ruido.)

Anclaje §A13 — matmul como búsqueda de embedding

Aunque este lab usa matrices aleatorias (N, N), recuerda el anclaje §A13 de theory/02-matmul-and-shapes.md: la búsqueda de embedding E @ one_hot(i) es un multiply matriz-vector con E.shape = (V, D) = (600, 64). A nuestra escala §A13 eso son 2 × 600 × 64 = 76,800 FLOPs por token — trivialmente dominado por el propio gather (E[i], cero FLOPs). A escala MiniGPT (V = 600, D = 64, B = 32, T = 16) por forward pass la proyección al vocabulario x @ E^T es 2 × 32 × 16 × 64 × 600 ≈ 40M FLOPs, que cae cómodamente en territorio BLAS. Saber dónde está tu operación en esta curva te dice si has de preocuparte por la sobrecarga de Python en absoluto.

Restricciones

  • Single-threaded. Pon OMP_NUM_THREADS=1, OPENBLAS_NUM_THREADS=1, MKL_NUM_THREADS=1 antes de lanzar bench.py. Registra las env vars en manifest.json. Los benchmarks multi-thread son Fase 35.
  • Governor de CPU performance + alimentación AC. Mismo protocolo que el lab 01 de la Fase 1.
  • Sin cython / numba / torch. NumPy puro / Python puro.
  • Mismo seed para todas las ejecuciones (para que los patrones de acceso a memoria no difieran entre métodos).

Escollos

  • np.matmul es no determinista entre ejecuciones por cantidades minúsculas (BLAS puede reordenar por motivos de cache). No esperes resultados bit-iguales entre ejecuciones; espera rtol=1e-5 entre métodos en la misma ejecución.
  • N=2048 podría OOM si asignas demasiados buffers. Reusa los mismos A, B, C entre iteraciones.
  • Naive en N=1024 tarda ~5-10 minutos por iteración en Python puro. Planifica según eso. Extrapolar desde N=512 es aceptable si lo documentas en el README.
  • np.einsum('ik,kj->ij', A, B) no es idéntico a A @ B en code path — ambos van a BLAS, pero por distintos caminos del optimizador. Ambas son mediciones válidas; registra ambas.
  • La primera llamada a una rutina BLAS en un proceso Python paga un coste de setup único. La iteración de warm-up lo gestiona; no midas arranque en frío.

Condiciones de parada

  • Las cuatro funciones implementadas en bench.py; la comprobación de acuerdo pasa con rtol=1e-5, atol=1e-6 en N=64.
  • El gráfico existe con las cuatro curvas + overlay del roofline.
  • El README responde a las cuatro preguntas de interpretación.
  • El speedup de np.matmul sobre matmul_naive en el N medido más grande es ≥ 50×. (Realísticamente será 10⁴-10⁵×; 50× es el piso.) Si no, depura.

Cuándo consultar solutions/

Tras commitear el gráfico + README. solutions/01-matmul-perf-ref.md (en la apertura de fase) recorre los números esperados en el i5-8250U y discute internals de BLAS a nivel curricular.

Siguiente lab: lab/02-svd-compression.md.