Algoritmo de multiplicación de matrices - Matrix multiplication algorithm

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 . Se han diseñado muchos algoritmos diferentes para multiplicar matrices en diferentes tipos de hardware, incluidos los 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 proporciona un algoritmo que toma tiempo del orden de $n 3$ operaciones de campo para multiplicar dos $n \times 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, 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 , 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 \times m$ y una matriz $B de$ $m$ $\times$ $p$ , entonces $C$ es una matriz de $n$ $\times$ $p$ con entradas

{\ Displaystyle c_ {ij} = \ sum _ {k = 1} ^ {m} a_ {ik} b_ {kj}.}

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 \leftarrow suma + A ik \times B kj$
  - Establecer $C ij \leftarrow suma$
Regresar $C$

Este algoritmo toma tiempo $Θ (nmp)$ (en notación asintótica ). 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 \times n$ , en cuyo caso el tiempo de ejecución es $Θ (n 3)$ , es decir, cúbico en el tamaño de la dimensión.

Comportamiento de la caché

Ilustración del orden principal de filas y columnas

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; qué orden es mejor 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, 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.

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 $\sqrt M$ por $\sqrt M$ :

Entrada: matrices $A$ y $B$
Sea $C$ una nueva matriz del tamaño apropiado
Elija un tamaño de mosaico $T = Θ (\sqrt 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 \leftarrow suma + A ik \times B kj$
        
        Establecer $C ij \leftarrow C ij + suma$
Regresar $C$

En el modelo de caché idealizado, este algoritmo solo incurre en $Θ (n 3 / b \sqrt M)$ fallas de caché; el divisor $b \sqrt 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é.

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.

{\ displaystyle C = {\ begin {pmatrix} C_ {11} & C_ {12} \\ C_ {21} & C_ {22} \\\ end {pmatrix}}, \, A = {\ begin {pmatrix} A_ { 11} & A_ {12} \\ A_ {21} & A_ {22} \\\ end {pmatrix}}, \, B = {\ begin {pmatrix} B_ {11} & B_ {12} \\ B_ {21} & B_ {22} \\\ end {pmatrix}},}

que funciona para todas las matrices cuadradas cuyas dimensiones son potencias de dos, es decir, las formas son $2 n \times 2 n$ para algunos $n$ . El producto de la matriz es ahora

{\ displaystyle {\ begin {pmatrix} C_ {11} & C_ {12} \\ C_ {21} & C_ {22} \\\ end {pmatrix}} = {\ begin {pmatrix} A_ {11} & A_ {12} \\ A_ {21} & A_ {22} \\\ end {pmatrix}} {\ begin {pmatrix} B_ {11} & B_ {12} \\ B_ {21} & B_ {22} \\\ end {pmatrix}} = {\ begin {pmatrix} A_ {11} B_ {11} + A_ {12} B_ {21} & A_ {11} B_ {12} + A_ {12} B_ {22} \\ A_ {21} B_ {11 } + A_ {22} B_ {21} & A_ {21} B_ {12} + A_ {22} B_ {22} \\\ end {pmatrix}}}

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

{\ Displaystyle T (1) = \ Theta (1);}

{\ Displaystyle T (n) = 8T (n / 2) + \ Theta (n ^ {2}),}

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.

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 divide las matrices en dos en lugar de cuatro submatrices, de la siguiente manera. 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 \times m$ , $B$ de tamaño $m \times 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:

{\ Displaystyle C = {\ begin {pmatrix} A_ {1} \\ A_ {2} \ end {pmatrix}} {B} = {\ begin {pmatrix} A_ {1} B \\ A_ {2} B \ end {pmatrix}}}

De lo contrario, si $max (n, m, p) = p$ , divida $B$ verticalmente:

{\ displaystyle C = A {\ begin {pmatrix} B_ {1} & B_ {2} \ end {pmatrix}} = {\ begin {pmatrix} AB_ {1} & AB_ {2} \ end {pmatrix}}}

De lo contrario, $max (n, m, p) = m$ . Divida $A$ verticalmente y $B$ horizontalmente:

{\ displaystyle C = {\ begin {pmatrix} A_ {1} & A_ {2} \ end {pmatrix}} {\ begin {pmatrix} B_ {1} \\ B_ {2} \ end {pmatrix}} = A_ { 1} B_ {1} + A_ {2} B_ {2}}

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é : no se requiere un 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é. (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 pérdidas de caché incurridas por este algoritmo, en una máquina con $M$ líneas de caché ideal, cada una de tamaño $b$ bytes, está limitada por

{\ Displaystyle \ Theta \ left (m + n + p + {\ frac {mn + np + mp} {b}} + {\ frac {mnp} {b {\ sqrt {M}}}} \ right)}

Algoritmos subcúbicos

Mejora de las estimaciones del exponente

ω a lo

largo del tiempo para la complejidad computacional de la multiplicación de matrices .

{\ Displaystyle O (n ^ {\ omega})}

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 \times 2$ que requiere solo 7 multiplicaciones (en lugar de las 8 habituales), a expensas de varias operaciones adicionales de suma y resta. Al aplicar esto de forma recursiva se obtiene 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, pero es más rápido en los casos en los que $n$ $> 100$ aproximadamente y aparece en varias bibliotecas, como BLAS . Es muy útil para matrices grandes sobre dominios exactos como campos finitos , donde la estabilidad numérica no es un problema. ${\ Displaystyle O (n ^ {\ log _ {2} 7}) \ approx O (n ^ {2.807})}$

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 multiplicación de matrices , generalmente denotado , es el número real más pequeño por el cual cualquier matriz sobre un campo se puede multiplicar mediante operaciones de campo. El mejor límite actual es el de Josh Alman y Virginia Vassilevska Williams . Este algoritmo, como todos los 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. La idea conceptual de estos algoritmos es similar al algoritmo de Strassen: una forma está diseñado para multiplicar dos $k$ $\times$ $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. ${\ Displaystyle \ omega}$ ${\ Displaystyle n \ times n}$ ${\ Displaystyle n ^ {\ omega + o (1)}}$ ${\ Displaystyle \ omega}$ ${\ Displaystyle \ omega <2.3728596}$

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

Paralelismo de memoria compartida

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

{\ Displaystyle {\ begin {pmatrix} A_ {11} B_ {11} + A_ {12} B_ {21} & A_ {11} B_ {12} + A_ {12} B_ {22} \\ A_ {21} B_ {11} + A_ {22} B_ {21} y A_ {21} B_ {12} + A_ {22} B_ {22} \\\ end {pmatrix}}}

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 :

Procedimiento $multiplicar (C, A, B)$ :

Caso base: si $n = 1$ , establezca $c 11 \leftarrow a 11 \times 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 \leftarrow 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 "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 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.

Multiplicación de matrices de bloques. En el algoritmo 2D, cada procesador es responsable de una submatriz de

C

. En el algoritmo 3D, cada par de submatrices de

A

y

B

que se multiplica se asigna a un procesador.

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 $\sqrt M / 3$ por $\sqrt M / 3$ , donde $M$ es el tamaño de la memoria rápida. El algoritmo ingenuo se usa luego sobre las matrices de bloques, calculando los productos de las submatrices completamente en la memoria rápida. Esto reduce el ancho de banda de comunicación a $O (n 3 / \sqrt M)$ , que es asintóticamente óptimo (para algoritmos que realizan cálculos $Ω (n 3)$ ).

En un entorno distribuido con $p$ procesadores dispuestos en una malla $\sqrt p$ por $\sqrt 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 / \sqrt p)$ palabras, lo cual es asintóticamente óptimo asumiendo que cada nodo almacena los elementos mínimos $O (n 2 / p)$ . 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. Este algoritmo transmite $O (n 2 / p 2/3)$ palabras por procesador, lo que es asintóticamente óptimo. 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. Los algoritmos "2.5D" proporcionan una compensación continua entre el uso de la memoria y el ancho de banda de comunicación. En entornos informáticos distribuidos modernos como MapReduce , se han desarrollado algoritmos de multiplicación especializados.

Algoritmos para mallas

Multiplicación de matrices completada en 2n-1 pasos para dos matrices n × n en una malla de cables cruzados.

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. La matriz estándar es ineficaz porque los datos de las dos matrices no llegan simultáneamente y debe 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. El rendimiento mejora aún más para cálculos repetidos que conducen a una eficiencia del 100%. La matriz de malla de alambres cruzados puede verse como un caso especial de una estructura de procesamiento no plana (es decir, multicapa).

Ver también

Referencias

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

Languages

In other projects

Algoritmo de multiplicación de matrices - Matrix multiplication algorithm

Contenido

Algoritmo iterativo

Comportamiento de la caché

Algoritmo de divide y vencerás

Matrices no cuadradas

Comportamiento de la caché

Algoritmos subcúbicos

Algoritmos paralelos y distribuidos

Paralelismo de memoria compartida

Algoritmos distribuidos y que evitan la comunicación

Algoritmos para mallas

Ver también

Referencias

Otras lecturas