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 las rutas a través de un gráfico . [1] Se han diseñado muchos algoritmos diferentes para multiplicar matrices en diferentes tipos de hardware, incluido el paralelo y el distribuido. sistemas, 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 proporciona un algoritmo que toma tiempo del orden de n 3 operaciones de campo para multiplicar dos n × n matrices sobre ese campo ( Θ ( n 3 ) en notación O grande ). Se conocen mejores límites asintóticos en el tiempo requerido para multiplicar matrices desde el algoritmo de Strassen en la década de 1960, pero aún se desconoce cuál es el tiempo óptimo (es decir, cuál es la complejidad computacional de la multiplicación de matrices ). A diciembre de 2020 [actualizar], el algoritmo de multiplicación de matrices con la mejor complejidad asintótica se ejecuta en el tiempo O ( n 2.3728596 ) , dado por Josh Alman y Virginia Vassilevska Williams , [2] [3] sin embargo, este algoritmo es un algoritmo galáctico debido a las grandes constantes y no se puede realizar en la práctica.
Algoritmo iterativo
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
De esto, un algoritmo simple se puede construir que los bucles a través de los índices i de 1 a n y j de 1 a p , el cálculo de la anterior usando un bucle anidado:
- Entrada: matrices A y B
- Sea C una nueva matriz del tamaño apropiado
- Para i de 1 an :
- Para j de 1 ap :
- Sea suma = 0
- Para k de 1 am :
- Establecer suma ← suma + A ik × B kj
- Establecer C ij ← suma
- Para j de 1 ap :
- Regresar C
Este algoritmo toma 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úbicos en el tamaño de la dimensión. . [4]
Comportamiento de la caché
Los tres bucles en la multiplicación iterativa de matrices se pueden intercambiar arbitrariamente entre sí sin que ello afecte a la corrección o al 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 caché del algoritmo; [1] cuál es 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 una caché totalmente asociativa que consta de M bytes yb bytes por línea de caché (es decir,METRO/Blíneas de caché), el algoritmo anterior es subó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 fallo de caché cuando se accede a un elemento de B . Esto significa que el algoritmo incurre en Θ ( n 3 ) fallas 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. [5]
La variante óptima del algoritmo iterativo para A y B en el diseño 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 : [5] [6]
- Entrada: matrices A y B
- Sea C una nueva matriz del tamaño apropiado
- Elija un tamaño de mosaico T = Θ ( √ M )
- Para I de 1 an en pasos de T :
- Para J de 1 ap en pasos de T :
- Para K de 1 am en pasos de T :
- Multiplica A I : I + T , K : K + T y B K : K + T , J : J + T en C I : I + T , J : J + T , es decir:
- Para i de I a min ( I + T , n ) :
- Para j de J a min ( J + T , p ) :
- Sea suma = 0
- Para k de K a min ( K + T , m ) :
- Establecer suma ← suma + A ik × B kj
- Establecer C ij ← C ij + suma
- Para j de J a min ( J + T , p ) :
- Para K de 1 am en pasos de T :
- Para J de 1 ap en pasos de T :
- Regresar C
En el modelo de caché idealizado, este algoritmo solo incurre en Θ ( n 3/b √ M) fallas 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 las fallas de caché. [5]
Algoritmo de divide y vencerás
Una alternativa al algoritmo iterativo es el algoritmo divide y vencerás para la multiplicación de matrices. Esto 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 algunos n . El producto de la matriz es ahora
que consta de ocho multiplicaciones de pares de submatrices, seguidas de un paso de suma. 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 su caso base.
La complejidad de este algoritmo en función de n viene dada por la recurrencia [4]
teniendo en cuenta las ocho llamadas recursivas en matrices de tamaño n / 2 y Θ ( n 2 ) para sumar los cuatro pares de matrices resultantes por elementos. La aplicación del teorema maestro para las recurrencias de divide y vencerás muestra que esta recursión tiene la solución Θ ( n 3 ) , lo mismo que el algoritmo iterativo. [4]
Matrices no cuadradas
Una variante de este algoritmo que funciona para matrices de formas arbitrarias y es más rápido en la práctica [5] divide las matrices en dos en lugar de cuatro submatrices, de la siguiente manera. [7] Dividir una matriz ahora significa dividirla en dos partes de igual tamaño, o lo más cerca posible de tamaños iguales en el caso de dimensiones impares.
- Entradas: matrices A de tamaño n × m , B de tamaño m × p .
- Caso base: si max ( n , m , p ) está por debajo de algún umbral, use una versión no enrollada del algoritmo iterativo.
- Casos recursivos:
- Si max ( n , m , p ) = n , divida A horizontalmente:
- De lo contrario, si max ( n , m , p ) = p , divida B verticalmente:
- De lo contrario, max ( n , m , p ) = m . Divida A verticalmente y B horizontalmente:
Comportamiento de la caché
La tasa de errores de caché de la multiplicación de matrices recursivas 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é : [7] 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 la caché son efectivamente dinámicos debido a que otros procesos ocupan espacio en la caché. [5] (El algoritmo iterativo simple también ignora la memoria caché, pero en la práctica es mucho más lento si el diseño de la matriz no se adapta al algoritmo).
El número de fallos de caché incurridos por este algoritmo, en una máquina con M líneas de caché ideal, cada uno de tamaño b bytes, está limitado por [7] : 13
Algoritmos subcúbicos
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 al que a menudo se hace referencia como "multiplicación rápida de matrices". Se basa en una forma de multiplicar dos matrices de 2 × 2 que requiere solo 7 multiplicaciones (en lugar de las 8 habituales), a expensas de varias operaciones adicionales de suma y resta. Aplicar esto de forma recursiva da un algoritmo con un costo multiplicativo de. El algoritmo de Strassen es más complejo y la estabilidad numérica se reduce en comparación con el algoritmo ingenuo, [8] pero es más rápido en los casos en los que n > 100 aproximadamente [1] y aparece en varias bibliotecas, como BLAS . [9] Es muy útil para matrices grandes sobre dominios exactos como campos finitos , donde la estabilidad numérica no es un problema.
Es una pregunta abierta en la informática teórica qué tan bien se puede mejorar el algoritmo de Strassen en términos de complejidad asintótica . El exponente de la multiplicación de matrices , generalmente denotado, es el número real más pequeño para el que cualquier La matriz sobre un campo se puede multiplicar usando operaciones de campo. El mejor destino actual es , por Josh Alman y Virginia Vassilevska Williams . [2] Este algoritmo, como todos los otros algoritmos recientes en esta línea de investigación, es una generalización del algoritmo Coppersmith-Winograd, que fue dado por Don Coppersmith y Shmuel Winograd en 1990. [10] La idea conceptual de estos algoritmos es similar al algoritmo de Strassen: se diseña 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. [11] [12]
Algoritmo Freivalds' es un simple algoritmo de Monte Carlo que, matrices dadas A , B y C , verifica en Θ ( n 2 ) de tiempo si AB = C .
Algoritmos paralelos y distribuidos
El algoritmo divide y vencerás esbozado anteriormente se puede paralelizar de dos formas para multiprocesadores de memoria compartida . Estos se basan en el hecho de que las ocho multiplicaciones de matrices recursivas en
se pueden realizar de forma independiente entre sí, al igual que las cuatro sumas (aunque el algoritmo necesita "unir" las multiplicaciones antes de hacer las sumas). Aprovechando el paralelismo completo del problema, se obtiene un algoritmo que se puede expresar en un pseudocódigo estilo fork-join : [13]
Procedimiento multiplicar ( C , A , B ) :
- Caso base: si n = 1 , establezca c 11 ← a 11 × b 11 (o multiplique una matriz de bloques pequeños).
- De lo contrario, asigne espacio para una nueva matriz T de forma n × n , luego:
- Divida A en A 11 , A 12 , A 21 , A 22 .
- Divida B en B 11 , B 12 , B 21 , B 22 .
- Divida C en C 11 , C 12 , C 21 , C 22 .
- Divida T en T 11 , T 12 , T 21 , T 22 .
- Ejecución paralela:
- Multiplicar por horquilla ( C 11 , A 11 , B 11 ) .
- Multiplicar por horquilla ( C 12 , A 11 , B 12 ) .
- Multiplicar por horquilla ( C 21 , A 21 , B 11 ) .
- Multiplicar por horquilla ( C 22 , A 21 , B 12 ) .
- Multiplicar por horquilla ( T 11 , A 12 , B 21 ) .
- Multiplicar por horquilla ( T 12 , A 12 , B 22 ) .
- Multiplicar por horquilla ( T 21 , A 22 , B 21 ) .
- Multiplicar por horquilla ( T 22 , A 22 , B 22 ) .
- Únase (espere a que se completen las bifurcaciones paralelas).
- agregar ( C , T ) .
- Deallocate T .
El procedimiento add ( C , T ) agrega T en C , elemento-sabio:
- Caso base: si n = 1 , establezca c 11 ← c 11 + t 11 (o haga un bucle corto, quizás desenrollado).
- De lo contrario:
- Divida C en C 11 , C 12 , C 21 , C 22 .
- Divida T en T 11 , T 12 , T 21 , T 22 .
- En paralelo:
- Horquilla añadir ( C 11 , T 11 ) .
- Horquilla añadir ( C 12 , T 12 ) .
- Horquilla añadir ( C 21 , T 21 ) .
- Horquilla añadir ( C 22 , T 22 ) .
- Únete .
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 previamente "bifurcados". 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 tanto 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 Θ ( n 2 ) , sin usar una matriz temporal. [13]
Algoritmos distribuidos y que evitan la comunicación
En 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 caché, mientras que en una máquina de múltiples nodos con memoria distribuida es la cantidad transferida entre los nodos; en cualquier caso, se denomina ancho de banda de comunicación . El algoritmo ingenuo que utiliza tres bucles anidados utiliza un ancho de banda de comunicación Ω ( n 3 ) .
El algoritmo de Cannon , también conocido como el algoritmo 2D , es un algoritmo que evita la comunicación 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. [14] El algoritmo ingenuo se usa luego sobre las matrices de bloques, calculando productos de submatrices completamente en 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 ) ). [15] [16]
En un entorno distribuido con p procesadores dispuestos en una malla √ p por √ p 2D, se puede asignar una submatriz del resultado a cada procesador, y el producto se puede calcular con cada procesador que transmite O ( n 2 / √ p ) palabras, lo cual es asintóticamente óptimo asumiendo que cada nodo almacena los elementos mínimos O ( n 2 / p ) . [16] Esto se puede mejorar con el algoritmo 3D, que organiza los procesadores en una malla de cubos 3D, asignando cada producto de dos submatrices de entrada a un solo procesador. Las submatrices de resultado se generan luego realizando una reducción en cada fila. [17] Este algoritmo transmite O ( n 2 / p 2/3 ) palabras por procesador, lo que es asintóticamente óptimo. [16] Sin embargo, esto requiere replicar cada elemento de la matriz de entrada p 1/3 veces, por lo que requiere un factor de p 1/3 más de 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. [17] Los algoritmos "2.5D" proporcionan una compensación continua entre el uso de memoria y el ancho de banda de comunicación. [18] En entornos informáticos distribuidos modernos como MapReduce , se han desarrollado algoritmos de multiplicación especializados. [19]
Algoritmos para mallas
Existe 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 de 2D Cannon , se 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. [20] 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 aún más rápido en una malla de alambre cruzado de dos capas, donde solo se necesitan 2 n -1 pasos. [21] El rendimiento mejora aún más para cálculos repetidos que conducen a una eficiencia del 100%. [22] La matriz de malla de cables cruzados puede verse como un caso especial de una estructura de procesamiento no plana (es decir, de varias capas). [23]
Ver también
- Complejidad computacional de operaciones matemáticas
- Complejidad computacional de la multiplicación de matrices
- Algoritmo CYK, algoritmo de §Valiant
- Multiplicación de cadenas de matrices
- Multiplicación dispersa matriz-vector
Referencias
- ↑ a b c d Skiena, Steven (2008). "Ordenar y buscar". El Manual de diseño de algoritmos . Saltador. págs. 45 –46, 401–403. doi : 10.1007 / 978-1-84800-070-4_4 . ISBN 978-1-84800-069-8.
- ^ a b Alman, Josh; Williams, Virginia Vassilevska (2020), "Un método láser refinado y una multiplicación de matrices más rápida", 32 ° Simposio anual ACM-SIAM sobre algoritmos discretos (SODA 2021) , arXiv : 2010.05846
- ^ Hartnett, Kevin. "La multiplicación de la matriz se acerca más a la meta mítica" . Revista Quanta . Consultado el 1 de abril de 2021 .
- ^ a b c Cormen, Thomas H .; Leiserson, Charles E .; Rivest, Ronald L .; Stein, Clifford (2009) [1990]. Introducción a los algoritmos (3ª ed.). MIT Press y McGraw-Hill. págs. 75–79. ISBN 0-262-03384-4.
- ^ a b c d e Amarasinghe, Saman; Leiserson, Charles (2010). "6.172 Ingeniería de Desempeño de Sistemas Software, Clase 8" . MIT OpenCourseWare . Instituto de Tecnología de Massachusetts . Consultado el 27 de enero de 2015 .
- ^ Lam, Monica S .; Rothberg, Edward E .; Wolf, Michael E. (1991). El rendimiento de la caché y las optimizaciones de los algoritmos bloqueados . Conf. Int. sobre soporte arquitectónico para lenguajes de programación y sistemas operativos (ASPLOS).
- ^ a b c Prokop, Harald (1999). Algoritmos ajenos a la caché (PDF) (Maestría). MIT.
- ^ Miller, Webb (1975), "Complejidad computacional y estabilidad numérica", SIAM News , 4 (2): 97–107, CiteSeerX 10.1.1.148.9947 , doi : 10.1137 / 0204009
- ^ Prensa, William H .; Flannery, Brian P .; Teukolsky, Saul A .; Vetterling, William T. (2007). Recetas numéricas: el arte de la informática científica (3ª ed.). Prensa de la Universidad de Cambridge . pag. 108 . ISBN 978-0-521-88068-8.
- ^ Calderero, Don; Winograd, Shmuel (1990), "Multiplicación de matrices mediante progresiones aritméticas" (PDF) , Journal of Symbolic Computation , 9 (3): 251, doi : 10.1016 / S0747-7171 (08) 80013-2
- ^ Iliopoulos, Costas S. (1989), "Límites de complejidad del peor caso en algoritmos para calcular la estructura canónica de grupos abelianos finitos y las formas normales de Hermite y Smith de una matriz entera" (PDF) , SIAM Journal on Computing , 18 (4 ): 658–669, CiteSeerX 10.1.1.531.9309 , doi : 10.1137 / 0218045 , MR 1004789 , archivado desde el original (PDF) en 2014-03-05 , obtenido 2015-01-16 ,
El algoritmo Coppersmith – Winograd no es práctico, debido a la gran constante oculta en el límite superior del número de multiplicaciones necesarias.
- ^ Robinson, Sara (noviembre de 2005), "Toward an Optimal Algorithm for Matrix Multiplication" (PDF) , SIAM News , 38 (9),
Incluso si alguien logra probar una de las conjeturas, demostrando así que ω = 2 - el producto de la corona Es poco probable que este enfoque 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 de tiempo sea evidente.
- ^ a b Randall, Keith H. (1998). Cilk: Computación multiproceso eficiente (PDF) (Ph.D.). Instituto de Tecnología de Massachusetts. págs. 54–57.
- ^ Lynn Elliot Cannon, Una computadora celular para implementar el algoritmo de filtro de Kalman , Informe técnico, Ph.D. Tesis, Universidad Estatal de Montana, 14 de julio de 1969.
- ^ Hong, JW; Kung, HT (1981). "Complejidad de E / S: el juego de guijarros rojo-azul" (PDF) . STOC '81: Actas del Decimotercer Simposio Anual de ACM sobre Teoría de la Computación : 326–333.
- ^ a b c Ironía, Dror; Toledo, Sivan; Tiskin, Alexander (septiembre de 2004). "Límites inferiores de comunicación para la multiplicación de matrices de memoria distribuida". J. Parallel Distrib. Computación . 64 (9): 1017–1026. CiteSeerX 10.1.1.20.7034 . doi : 10.1016 / j.jpdc.2004.03.021 .
- ^ a b Agarwal, RC; Balle, SM; Gustavson, FG; Joshi, M .; Palkar, P. (septiembre de 1995). "Un enfoque tridimensional para la multiplicación de matrices paralelas". IBM J. Res. Dev . 39 (5): 575–582. CiteSeerX 10.1.1.44.3404 . doi : 10.1147 / rd.395.0575 .
- ^ Solomonik, Edgar; Demmel, James (2011). "Algoritmos de factorización y multiplicación de matrices 2.5D paralelas óptimas de comunicación". Actas de la 17ª Conferencia Internacional sobre Procesamiento Paralelo . Parte II: 90–109.
- ^ Bosagh Zadeh, Reza; Carlsson, Gunnar (2013). "Cuadrado de matriz independiente de dimensión utilizando MapReduce" (PDF) . arXiv : 1304.1467 . Código Bibliográfico : 2013arXiv1304.1467B . Consultado el 12 de julio de 2014 . Cite journal requiere
|journal=
( ayuda ) - ^ Bae, SE; Shinn, T.-W .; Takaoka, T. (2014). "Un algoritmo paralelo más rápido para la multiplicación de matrices en una matriz de malla" . Procedia Informática . 29 : 2230–2240. doi : 10.1016 / j.procs.2014.05.208 .
- ^ Kak, S (1988). "Una matriz de malla de dos capas para la multiplicación de matrices". Computación paralela . 6 (3): 383–385. CiteSeerX 10.1.1.88.8527 . doi : 10.1016 / 0167-8191 (88) 90078-6 .
- ^ Kak, S. (2014) Eficiencia de la multiplicación de matrices en la matriz de malla de alambre cruzado. https://arxiv.org/abs/1411.3273
- ^ Kak, S (1988). "Computación de matriz multicapa". Ciencias de la información . 45 (3): 347–365. CiteSeerX 10.1.1.90.4753 . doi : 10.1016 / 0020-0255 (88) 90010-2 .
Otras lecturas
- Buttari, Alfredo; Langou, Julien; Kurzak, Jakub; Dongarra, Jack (2009). "Una clase de algoritmos de álgebra lineal en mosaico paralelo para arquitecturas multinúcleo". Computación paralela . 35 : 38–53. arXiv : 0709.1272 . doi : 10.1016 / j.parco.2008.10.002 . S2CID 955 .
- Goto, Kazushige; van de Geijn, Robert A. (2008). "Anatomía de la multiplicación de matrices de alto rendimiento". Transacciones ACM en software matemático . 34 (3): 1–25. CiteSeerX 10.1.1.140.3583 . doi : 10.1145 / 1356052.1356053 . S2CID 9359223 .
- Van Zee, Field G .; van de Geijn, Robert A. (2015). "BLIS: un marco para instanciar rápidamente la funcionalidad BLAS". Transacciones ACM en software matemático . 41 (3): 1–33. doi : 10.1145 / 2764454 . S2CID 1242360 .
- Cómo optimizar GEMM