Debido a que la multiplicación de matrices es una operación tan central en muchos algoritmos numéricos , se ha invertido mucho trabajo en hacer que los algoritmos de multiplicación de matrices sean eficientes. Las aplicaciones de la multiplicación de matrices en problemas computacionales se encuentran en muchos campos, incluida la computación científica y el reconocimiento de patrones y en problemas aparentemente no relacionados, como contar los caminos a través de un gráfico . [1] Se han diseñado muchos algoritmos diferentes para multiplicar matrices en diferentes tipos de hardware, incluidos sistemas paralelos y distribuidos , donde el trabajo computacional se distribuye en múltiples procesadores (quizás en una red).
La aplicación directa de la definición matemática de multiplicación de matrices da como resultado un algoritmo que tarda un tiempo del orden de n 3 operaciones de campo en multiplicar dos matrices n × n sobre ese campo ( Θ( n 3 ) en notación O grande ). Se conocen mejores límites asintóticos del tiempo necesario para multiplicar matrices desde el algoritmo de Strassen en la década de 1960, pero el tiempo óptimo (es decir, la complejidad computacional de la multiplicación de matrices ) sigue siendo desconocido. A partir de abril de 2024 [actualizar], el mejor límite anunciado para la complejidad asintótica de un algoritmo de multiplicación de matrices es el tiempo O( n 2,371552 ) , dado por Williams , Xu, Xu y Zhou. [2] [3] Esto mejora el límite de tiempo O( n 2,3728596 ) , dado por Alman y Williams. [4] [5] Sin embargo, este algoritmo es un algoritmo galáctico debido a las grandes constantes y no se puede realizar de manera práctica.
La definición de multiplicación de matrices es que si C = AB para una matriz A de n × m y una matriz B de m × p , entonces C es una matriz de n × p con entradas
A partir de esto, se puede construir un algoritmo simple que recorre los índices i desde 1 hasta n y j desde 1 hasta p , calculando lo anterior mediante un bucle anidado:
Este algoritmo toma un tiempo Θ( nmp ) (en notación asintótica ). [1] Una simplificación común para el propósito del análisis de algoritmos es asumir que las entradas son todas matrices cuadradas de tamaño n × n , en cuyo caso el tiempo de ejecución es Θ( n 3 ) , es decir, cúbico en el tamaño de la dimensión. [6]
Los tres bucles en la multiplicación iterativa de matrices se pueden intercambiar arbitrariamente entre sí sin que esto afecte la corrección o el tiempo de ejecución asintótico. Sin embargo, el orden puede tener un impacto considerable en el rendimiento práctico debido a los patrones de acceso a la memoria y al uso de la memoria caché del algoritmo; [1] el mejor orden también depende de si las matrices se almacenan en orden de fila principal, orden de columna principal o una combinación de ambos.
En particular, en el caso idealizado de un caché completamente asociativo que consta de M bytes y b bytes por línea de caché (es decir ,METRO/b líneas de caché), el algoritmo anterior no es óptimo para A y B almacenados en orden de fila principal. Cuando n > METRO/b , cada iteración del bucle interno (un barrido simultáneo a través de una fila de A y una columna de B ) incurre en un error de caché al acceder a un elemento de B . Esto significa que el algoritmo incurre en Θ( n 3 ) errores de caché en el peor de los casos. A partir de 2010[actualizar], la velocidad de las memorias en comparación con la de los procesadores es tal que los errores de caché, en lugar de los cálculos reales, dominan el tiempo de ejecución para matrices considerables. [7]
La variante óptima del algoritmo iterativo para A y B en disposición de filas principales es una versión en mosaico , donde la matriz se divide implícitamente en mosaicos cuadrados de tamaño √ M por √ M : [7] [8]
En el modelo de caché idealizado, este algoritmo solo incurre en Θ( número 3/b√M ) errores de caché; el divisor b √ M equivale a varios órdenes de magnitud en las máquinas modernas , de modo que los cálculos reales dominan el tiempo de ejecución, en lugar de los errores de caché. [7]
Una alternativa al algoritmo iterativo es el algoritmo de divide y vencerás para la multiplicación de matrices. Este se basa en la partición de bloques.
que funciona para todas las matrices cuadradas cuyas dimensiones son potencias de dos, es decir, las formas son 2 n × 2 n para algún n . El producto matricial es ahora
que consiste en ocho multiplicaciones de pares de submatrices, seguidas de un paso de adición. El algoritmo divide y vencerás calcula las multiplicaciones más pequeñas de forma recursiva , utilizando la multiplicación escalar c 11 = a 11 b 11 como caso base.
La complejidad de este algoritmo en función de n viene dada por la recurrencia [6]
teniendo en cuenta las ocho llamadas recursivas a matrices de tamaño n /2 y Θ( n 2 ) para sumar los cuatro pares de matrices resultantes elemento por elemento. La aplicación del teorema maestro para recurrencias de divide y vencerás muestra que esta recursión tiene la solución Θ( n 3 ) , la misma que el algoritmo iterativo. [6]
Una variante de este algoritmo que funciona para matrices de formas arbitrarias y es más rápida en la práctica [7] divide las matrices en dos en lugar de cuatro submatrices, de la siguiente manera. [9] Dividir una matriz ahora significa dividirla en dos partes de igual tamaño, o lo más cercanas a tamaños iguales como sea posible en el caso de dimensiones impares.
La tasa de errores de caché de la multiplicación de matrices recursiva es la misma que la de una versión iterativa en mosaico , pero a diferencia de ese algoritmo, el algoritmo recursivo no tiene en cuenta la caché : [9] no se requiere ningún parámetro de ajuste para obtener un rendimiento óptimo de la caché, y se comporta bien en un entorno de multiprogramación donde los tamaños de caché son efectivamente dinámicos debido a que otros procesos ocupan espacio de caché. [7] (El algoritmo iterativo simple también no tiene en cuenta la caché, pero es mucho más lento en la práctica si el diseño de la matriz no está adaptado al algoritmo).
La cantidad de errores de caché que se producen con este algoritmo, en una máquina con M líneas de caché ideal, cada una de tamaño b bytes, está limitada por [9] : 13
Existen algoritmos que proporcionan mejores tiempos de ejecución que los sencillos. El primero en ser descubierto fue el algoritmo de Strassen , ideado por Volker Strassen en 1969 y a menudo denominado "multiplicación rápida de matrices". Se basa en una forma de multiplicar dos matrices de 2 × 2 que requieren solo 7 multiplicaciones (en lugar de las 8 habituales), a expensas de varias operaciones de suma y resta adicionales. Al aplicar esto de forma recursiva, se obtiene un algoritmo con un coste multiplicativo de . El algoritmo de Strassen es más complejo y la estabilidad numérica se reduce en comparación con el algoritmo ingenuo, [10] pero es más rápido en los casos en los que n > 100 o así [1] y aparece en varias bibliotecas, como BLAS . [11] Es muy útil para matrices grandes sobre dominios exactos como los campos finitos , donde la estabilidad numérica no es un problema.
Dado que el algoritmo de Strassen se utiliza actualmente en software numérico práctico y sistemas de álgebra computacional, mejorar las constantes ocultas en la notación O grande tiene sus méritos. A continuación se incluye una tabla que compara aspectos clave de la versión mejorada basada en la multiplicación recursiva de matrices de bloques 2x2 mediante multiplicaciones de matrices de 7 bloques. Como es habitual, se indican las dimensiones de la matriz y se designa el tamaño de la memoria.
Se sabe que un algoritmo similar al de Strassen con un paso de matriz de bloques 2x2 requiere al menos 7 multiplicaciones de matrices de bloques. En 1976, Probert [16] demostró que un algoritmo de este tipo requiere al menos 15 adiciones (incluidas las restas), sin embargo, una suposición oculta era que los bloques y la matriz de bloques 2x2 están representados en la misma base. Karstadt y Schwartz calcularon en diferentes bases y cambiaron 3 adiciones por transformaciones de base menos costosas. También demostraron que no se puede ir por debajo de 12 adiciones por paso utilizando diferentes bases. En un trabajo posterior, Beniamini et al. [17] aplicaron este truco de cambio de base a descomposiciones más generales que las matrices de bloques 2x2 y mejoraron la constante principal para sus tiempos de ejecución.
Es una pregunta abierta en la informática teórica hasta qué punto se puede mejorar el algoritmo de Strassen en términos de complejidad asintótica . El exponente de multiplicación de matrices , normalmente denotado como , es el número real más pequeño para el que cualquier matriz sobre un cuerpo se puede multiplicar entre sí mediante operaciones de cuerpo. La mejor cota actual de es , por Williams , Xu, Xu y Zhou. [2] [4] Este algoritmo, como todos los demás algoritmos recientes en esta línea de investigación, es una generalización del algoritmo Coppersmith-Winograd, que fue propuesto por Don Coppersmith y Shmuel Winograd en 1990. [18] La idea conceptual de estos algoritmos es similar al algoritmo de Strassen: se idea una forma de multiplicar dos k × k -matrices con menos de k 3 multiplicaciones, y esta técnica se aplica de forma recursiva. Sin embargo, el coeficiente constante oculto por la notación Big O es tan grande que estos algoritmos solo valen la pena para matrices que son demasiado grandes para manejarlas en las computadoras actuales. [19] [20]
El algoritmo de Freivalds es un algoritmo de Monte Carlo simple que , dadas las matrices A , B y C , verifica en tiempo Θ( n 2 ) si AB = C.
En 2022, DeepMind presentó AlphaTensor, una red neuronal que utilizó una analogía de juego de un solo jugador para inventar miles de algoritmos de multiplicación de matrices, incluidos algunos descubiertos previamente por humanos y otros que no. [21] Las operaciones se restringieron al campo fundamental no conmutativo (aritmética normal) y al campo finito (aritmética mod 2). El mejor algoritmo "práctico" (descomposición explícita de bajo rango de un tensor de multiplicación de matrices) encontrado funcionó en O(n 2.778 ). [22] Encontrar descomposiciones de bajo rango de tales tensores (y más allá) es NP-hard; la multiplicación óptima incluso para matrices 3x3 sigue siendo desconocida , incluso en el campo conmutativo. [22] En matrices 4x4, AlphaTensor descubrió inesperadamente una solución con 47 pasos de multiplicación, una mejora con respecto a los 49 requeridos con el algoritmo de Strassen de 1969, aunque restringida a la aritmética mod 2. De manera similar, AlphaTensor resolvió matrices de 5x5 con 96 pasos en lugar de los 98 de Strassen. Basándose en el sorprendente descubrimiento de que tales mejoras existen, otros investigadores pudieron encontrar rápidamente un algoritmo 4x4 independiente similar y ajustaron por separado el algoritmo 5x5 de 96 pasos de Deepmind a 95 pasos en aritmética módulo 2 y a 97 [23] en aritmética normal. [24] Algunos algoritmos nunca se habían descubierto antes, por ejemplo, (4, 5, 5) se mejoró de 80 a 76 en aritmética normal y módulo 2.
El algoritmo de dividir y vencer esbozado anteriormente se puede paralelizar de dos maneras para multiprocesadores de memoria compartida . Estas se basan en el hecho de que las ocho multiplicaciones de matrices recursivas en
se pueden realizar independientemente una de la otra, al igual que las cuatro sumas (aunque el algoritmo necesita "unir" las multiplicaciones antes de realizar las sumas). Al explotar el paralelismo completo del problema, se obtiene un algoritmo que se puede expresar en pseudocódigo de estilo fork-join : [25]
Procedimiento multiplicar( C , A , B ) :
El procedimiento add( C , T ) agrega T a C , elemento por elemento:
Aquí, fork es una palabra clave que indica que un cálculo puede ejecutarse en paralelo con el resto de la llamada a la función, mientras que join espera a que se completen todos los cálculos "bifurcados" previamente. La partición logra su objetivo únicamente mediante la manipulación del puntero.
Este algoritmo tiene una longitud de ruta crítica de Θ(log 2 n ) pasos, lo que significa que lleva ese tiempo en una máquina ideal con un número infinito de procesadores; por lo tanto, tiene una aceleración máxima posible de Θ( n 3 /log 2 n ) en cualquier computadora real. El algoritmo no es práctico debido al costo de comunicación inherente al movimiento de datos hacia y desde la matriz temporal T , pero una variante más práctica logra una aceleración de Θ( n 2 ) , sin usar una matriz temporal. [25]
En las arquitecturas modernas con memoria jerárquica, el costo de cargar y almacenar elementos de la matriz de entrada tiende a dominar el costo de la aritmética. En una sola máquina, esta es la cantidad de datos transferidos entre la RAM y la memoria caché, mientras que en una máquina multinodo con memoria distribuida es la cantidad transferida entre nodos; en ambos casos se denomina ancho de banda de comunicación . El algoritmo ingenuo que utiliza tres bucles anidados utiliza un ancho de banda de comunicación de Ω( n 3 ) .
El algoritmo de Cannon , también conocido como algoritmo 2D , es un algoritmo que evita la comunicación y que divide cada matriz de entrada en una matriz de bloques cuyos elementos son submatrices de tamaño √ M /3 por √ M /3 , donde M es el tamaño de la memoria rápida. [26] Luego, se utiliza el algoritmo ingenuo sobre las matrices de bloques, calculando productos de submatrices completamente en la memoria rápida. Esto reduce el ancho de banda de comunicación a O ( n 3 / √ M ) , que es asintóticamente óptimo (para algoritmos que realizan cálculos Ω( n 3 ) ). [27] [28]
En un entorno distribuido con p procesadores dispuestos en una malla 2D de √ p por √ p , se puede asignar una submatriz del resultado a cada procesador, y el producto se puede calcular con cada procesador transmitiendo O ( n 2 / √ p ) palabras, lo que es asintóticamente óptimo suponiendo que cada nodo almacena los O ( n 2 / p ) elementos mínimos. [28] Esto se puede mejorar con el algoritmo 3D, que organiza los procesadores en una malla de cubo 3D, asignando cada producto de dos submatrices de entrada a un solo procesador. Las submatrices de resultado se generan luego realizando una reducción sobre cada fila. [29] Este algoritmo transmite O ( n 2 / p 2/3 ) palabras por procesador, lo que es asintóticamente óptimo. [28] Sin embargo, esto requiere replicar cada elemento de la matriz de entrada p 1/3 veces, y por lo tanto requiere un factor de p 1/3 más memoria de la necesaria para almacenar las entradas. Este algoritmo se puede combinar con Strassen para reducir aún más el tiempo de ejecución. [29] Los algoritmos "2.5D" proporcionan un equilibrio continuo entre el uso de la memoria y el ancho de banda de comunicación. [30] En entornos informáticos distribuidos modernos como MapReduce , se han desarrollado algoritmos de multiplicación especializados. [31]
Hay una variedad de algoritmos para la multiplicación en mallas . Para la multiplicación de dos n × n en una malla bidimensional estándar utilizando el algoritmo 2D de Cannon , uno puede completar la multiplicación en 3 n -2 pasos aunque esto se reduce a la mitad de este número para cálculos repetidos. [32] La matriz estándar es ineficiente porque los datos de las dos matrices no llegan simultáneamente y deben rellenarse con ceros.
El resultado es incluso más rápido en una malla de alambre cruzado de dos capas, donde solo se necesitan 2 n -1 pasos. [33] El rendimiento mejora aún más para los cálculos repetidos, lo que lleva a una eficiencia del 100%. [34] La matriz de malla de alambre cruzado puede verse como un caso especial de una estructura de procesamiento no plana (es decir, multicapa). [35]
El algoritmo Coppersmith–Winograd no es práctico, debido a la constante oculta muy grande en el límite superior del número de multiplicaciones requeridas.
Incluso si alguien logra demostrar una de las conjeturas (demostrando así que
ω
= 2
), es poco probable que el enfoque del producto de corona sea aplicable a los grandes problemas de matrices que surgen en la práctica. [...] las matrices de entrada deben ser astronómicamente grandes para que la diferencia en el tiempo sea evidente.