Introducción

La multiplicación de matrices es una operación fundamental en áreas como la ciencia de datos, visión por computadora, simulaciones físicas y machine learning. En Python existen múltiples formas de llevarla a cabo, desde implementaciones puras en bucles for hasta librerías altamente optimizadas que aprovechan SIMD, multihilos y GPUs.

Conceptos Matemáticos Básicos

Dados dos matrices A (m × n) y B (n × p), su producto C = A·B es una matriz m × p cuyas entradas se calculan como:

c_{ij} = \sum_{k=1}^{n} a_{ik} \times b_{kj}

Este cálculo tiene una complejidad teórica O(m·n·p). Cuando m = n = p = N, el coste es O(N³), aunque existen algoritmos que reducen el exponente.

Algoritmos de Multiplicación

def mat_mul_naive(A, B):
    m, n = len(A), len(A[0])
    n2, p = len(B), len(B[0])
    assert n == n2, "Dimensiones incompatibles"
    C = [[0] * p for _ in range(m)]
    for i in range(m):
        for j in range(p):
            suma = 0
            for k in range(n):
                suma += A[i][k] * B[k][j]
            C[i][j] = suma
    return C

import numpy as np
A = np.random.rand(500, 500)
B = np.random.rand(500, 500)
C = A @ B   # o np.dot(A, B)

Algoritmo de Strassen (Dividir y Conquistar)

Strassen reduce la complejidad a O(N^{2.81}) al usar 7 multiplicaciones recursivas en lugar de 8. Es útil para tamaños de potencia de 2 y cuando se busca reducir la carga de cálculo en CPU.

def strassen(A, B):
    """Multiplicación de matrices usando el algoritmo de Strassen.
    Requiere que A y B sean cuadradas y de dimensión potencia de 2.
    """
    import numpy as np
    n = A.shape[0]
    if n <= 64:  # umbral: usar NumPy para sub‑matrices pequeñas
        return A @ B
    # Dividir en sub‑bloques
    mid = n // 2
    a11, a12 = A[:mid, :mid], A[:mid, mid:]
    a21, a22 = A[mid:, :mid], A[mid:, mid:]
    b11, b12 = B[:mid, :mid], B[:mid, mid:]
    b21, b22 = B[mid:, :mid], B[mid:, mid:]
    # 7 productos recursivos
    m1 = strassen(a11 + a22, b11 + b22)
    m2 = strassen(a21 + a22, b11)
    m3 = strassen(a11, b12 - b22)
    m4 = strassen(a22, b21 - b11)
    m5 = strassen(a11 + a12, b22)
    m6 = strassen(a21 - a11, b11 + b12)
    m7 = strassen(a12 - a22, b21 + b22)
    # Reconstruir C
    c11 = m1 + m4 - m5 + m7
    c12 = m3 + m5
    c21 = m2 + m4
    c22 = m1 - m2 + m3 + m6
    C = np.vstack((np.hstack((c11, c12)), np.hstack((c21, c22))))
    return C

Consideraciones de rendimiento:

• El algoritmo introduce sobrecarga de memoria (creación de sub‑matrices).
• Para N < 256 suele ser más lento que BLAS.
• Escala mejor en entornos donde el número de núcleos es limitado pero la latencia de memoria es alta.

Multiplicación en GPU con CuPy

Cuando se dispone de una GPU compatible con CUDA, CuPy ofrece una API idéntica a NumPy pero ejecuta los cálculos en la tarjeta gráfica.

import cupy as cp
A = cp.random.rand(2000, 2000, dtype=cp.float32)
B = cp.random.rand(2000, 2000, dtype=cp.float32)
C = A @ B   # Operación totalmente en GPU
# Transferir a CPU si es necesario
C_cpu = cp.asnumpy(C)

Ventajas: órdenes de magnitud más rápido para N > 1000.
Desventajas: requiere hardware CUDA y la gestión de la memoria GPU.

Benchmark Comparativo (CPU vs GPU)

El benchmark muestra que, a partir de N≈512, la solución basada en BLAS ya supera al algoritmo puro y que la GPU ofrece mejoras sustanciales para tamaños mayores.

Buenas Prácticas y Optimización

• Usar siempre NumPy o librerías basadas en BLAS cuando la precisión de punto flotante es suficiente.
• Evitar bucles anidados en Python puro; si es necesario, emplear numba o Cython para compilar.
• Prefijar dimensiones a potencias de 2 si se desea aplicar Strassen.
• Monitorizar el uso de memoria con tracemalloc o herramientas del sistema; la multiplicación de matrices grandes puede consumir GB de RAM.
• Utilizar tipos de datos adecuados: float32 en GPUs para reducir ancho de banda, float64 en CPU si se necesita precisión doble.

Resolución de Problemas Comunes (Troubleshooting)

1. Dimensiones incompatibles: verifica que el número de columnas de A coincida con el número de filas de B. Usa A.shape y B.shape para depuración.
2. MemoryError en Python puro: divide la operación en bloques ("block matrix multiplication") o cambia a NumPy.
3. Rendimiento inesperadamente bajo con NumPy: asegúrate de que la librería subyacente (OpenBLAS, MKL) está habilitada y que el número de hilos está configurado (por ejemplo, export OMP_NUM_THREADS=8).
4. Errores de CUDA con CuPy: verifica la versión del driver, la compatibilidad de la tarjeta y la instalación de cupy-cuda12x correspondiente.

Conclusiones

La multiplicación de matrices es un caso clásico donde la elección del algoritmo y la herramienta determinan la diferencia entre minutos y horas de cálculo. En la práctica, la regla de oro es:

"Empieza con numpy.dot o el operador @. Solo recurre a implementaciones personalizadas (Strassen, Numba, CuPy) cuando la métrica de rendimiento lo justifique y tengas los recursos adecuados."

Con los ejemplos y comparativas presentados, ahora puedes seleccionar la solución óptima para tu caso de uso, ya sea una aplicación de ciencia de datos que procesa cientos de miles de filas o un motor de IA que requiere cálculos masivos en GPU.