Algoritmo de multiplicação de matriz - Matrix multiplication algorithm

Como a multiplicação de matrizes é uma operação central em muitos algoritmos numéricos , muito trabalho foi investido para tornar os algoritmos de multiplicação de matrizes eficientes. Aplicações de multiplicação de matrizes em problemas computacionais são encontradas em muitos campos, incluindo computação científica e reconhecimento de padrões e em problemas aparentemente não relacionados, como contagem de caminhos em um gráfico . Muitos algoritmos diferentes foram projetados para multiplicar matrizes em diferentes tipos de hardware, incluindo sistemas paralelos e distribuídos , onde o trabalho computacional é distribuído por vários processadores (talvez em uma rede).

Aplicando diretamente a definição matemática de multiplicação de matrizes dá um algoritmo que leva tempo na ordem de $n 3$ campo operações para multiplicar dois $n \times n$ matrizes sobre esse campo ( $Θ (n 3)$ na notação O grande ). Melhores limites assintóticos no tempo necessário para multiplicar matrizes são conhecidos desde o algoritmo de Strassen na década de 1960, mas ainda não se sabe qual é o tempo ótimo (isto é, qual é a complexidade computacional da multiplicação de matrizes ). Em dezembro de 2020, o algoritmo de multiplicação de matrizes com melhor complexidade assintótica é executado no tempo $O (n 2,3728596)$ , fornecido por Josh Alman e Virginia Vassilevska Williams , no entanto, este algoritmo é um algoritmo galáctico devido às grandes constantes e não pode ser realizado na prática.

Algoritmo iterativo

A definição de multiplicação de matrizes é que se $C = AB$ para uma matriz $A$ $n \times m$ e uma matriz $B$ $m$ $\times$ $p$ , então $C$ é uma matriz $n$ $\times$ $p$ com entradas

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

A partir disso, um algoritmo simples pode ser construído que faz um loop sobre os índices $i$ de 1 a $n$ e $j$ de 1 a $p$ , computando o acima usando um loop aninhado:

Entrada: matrizes $A$ e $B$
Seja $C$ uma nova matriz de tamanho apropriado
Para i de 1 a n :
- Para j de 1 a p :
  - Deixe $somar = 0$
  - Para k de 1 a m :
    - Definir $soma \leftarrow soma + A ik \times B kj$
  - Definir $C ij \leftarrow soma$
Return $C$

Este algoritmo leva tempo $Θ (nmp)$ (em notação assintótica ). Uma simplificação comum para o propósito de análise de algoritmos é assumir que as entradas são todas matrizes quadradas de tamanho $n \times n$ , caso em que o tempo de execução é $Θ (n 3)$ , ou seja, cúbico no tamanho da dimensão.

Comportamento do cache

Ilustração da ordem principal de linha e coluna

Os três loops na multiplicação da matriz iterativa podem ser trocados arbitrariamente entre si sem um efeito na exatidão ou no tempo de execução assintótico. No entanto, o pedido pode ter um impacto considerável no desempenho prático devido aos padrões de acesso à memória e ao uso de cache do algoritmo; qual ordem é a melhor também depende se as matrizes são armazenadas na ordem principal da linha, ordem principal da coluna ou uma combinação de ambas.

Em particular, no caso idealizado de um cache totalmente associativo consistindo em $M$ bytes $eb$ bytes por linha de cache (ou seja,M/blinhas de cache), o algoritmo acima é abaixo do ideal para $A$ e $B$ armazenados na ordem da linha principal. Quando $n > M / b$ , Cada iteração do circuito interno (uma varredura simultânea por meio de uma linha de $um$ e uma coluna de $B$ ) incorre numa cache miss o acesso a um elemento de $B$ . Isso significa que o algoritmo incorre em $Θ (n 3)$ falhas de cache no pior caso. A partir de 2010, a velocidade das memórias em comparação com a dos processadores é tal que as perdas de cache, ao invés dos cálculos reais, dominam o tempo de execução para matrizes consideráveis.

A variante ótima do algoritmo iterativo para $A$ e $B$ no layout de linha maior é uma versão lado a lado , onde a matriz é implicitamente dividida em blocos quadrados de tamanho $\sqrt M$ por $\sqrt M$ :

Entrada: matrizes $A$ e $B$
Seja $C$ uma nova matriz de tamanho apropriado
Escolha um tamanho de ladrilho $T = Θ (\sqrt M)$
Para I de 1 a n nas etapas de T :
- Para J de 1 a p nas etapas de T :
  - Para K de 1 a m em etapas de T :
    - Multiplique $A I : I + T, K : K + T$ e $B K : K + T, J : J + T$ em $C I : I + T, J : J + T$ , isto é:
    - Para i de I a min ( I + T , n ) :
      - Para $j$ de $J$ a $min (J + T, p)$ :
        
        Deixe $somar = 0$
        
        Para $k$ de $K$ a $min (K + T, m)$ :
        Definir $soma \leftarrow soma + A ik \times B kj$
        
        Definir $C ij \leftarrow C ij + soma$
Return $C$

No modelo de cache idealizado, este algoritmo incorre apenas em $Θ (n 3 / b \sqrt M)$ falhas de cache; o divisor $b \sqrt M$ equivale a várias ordens de magnitude nas máquinas modernas, de modo que os cálculos reais dominam o tempo de execução, em vez de erros do cache.

Algoritmo de divisão e conquista

Uma alternativa ao algoritmo iterativo é o algoritmo de divisão e conquista para multiplicação de matrizes. Isso depende do particionamento de blocos

{\ 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 {pmatriz}}, \, B = {\ begin {pmatrix} B_ {11} & B_ {12} \\ B_ {21} & B_ {22} \\\ fim {pmatriz}},}

que funciona para todas as matrizes quadradas cujas dimensões são potências de dois, ou seja, as formas são $2 n \times 2 n$ para alguns $n$ . O produto da matriz é agora

{\ displaystyle {\ begin {pmatrix} C_ {11} & C_ {12} \\ C_ {21} & C_ {22} \\\ end {pmatrix}} = {\ begin {pmatrix} A_ {11} & A_ {12} \\ A_ {21} & A_ {22} \\\ end {pmatriz}} {\ 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} e A_ {21} B_ {12} + A_ {22} B_ {22} \\\ end {pmatriz}}}

que consiste em oito multiplicações de pares de submatrizes, seguidas por uma etapa de adição. O algoritmo de divisão e conquista calcula as multiplicações menores recursivamente , usando a multiplicação escalar $c 11 = a 11 b 11$ como seu caso base.

A complexidade deste algoritmo em função de $n$ é dada pela recorrência

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

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

levando em consideração as oito chamadas recursivas em matrizes de tamanho $n / 2$ e $Θ (n 2)$ para somar os quatro pares de matrizes resultantes elemento a elemento. A aplicação do teorema mestre para recorrências de divisão e conquista mostra que essa recursão tem a solução $Θ (n 3)$ , a mesma do algoritmo iterativo.

Matrizes não quadradas

Uma variante desse algoritmo, que funciona para matrizes de formas arbitrárias e é mais rápida na prática, divide as matrizes em duas em vez de quatro submatrizes, como segue. Dividir uma matriz agora significa dividi-la em duas partes de tamanho igual ou o mais próximo possível de tamanhos iguais no caso de dimensões ímpares.

Entradas: matrizes $A$ de tamanho $n \times m$ , $B$ de tamanho $m \times p$ .
Caso base: se $max (n, m, p)$ estiver abaixo de algum limite, use uma versão não rolada do algoritmo iterativo.
Casos recursivos:

Se $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 \ fim {pmatrix}}}

Caso contrário, se $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}}}

Caso contrário, $max (n, m, p) = m$ . Divida $A$ verticalmente e $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}}

Comportamento do cache

A taxa de perda de cache de multiplicação de matriz recursiva é a mesma que a de uma versão iterativa lado a lado , mas ao contrário desse algoritmo, o algoritmo recursivo é inconsciente do cache : não há parâmetro de ajuste necessário para obter o desempenho ideal do cache e se comporta bem em um ambiente de multiprogramação onde os tamanhos de cache são efetivamente dinâmicos devido a outros processos que ocupam espaço de cache. (O algoritmo iterativo simples também ignora o cache, mas é muito mais lento na prática se o layout da matriz não estiver adaptado ao algoritmo.)

O número de falhas de cache incorridas por este algoritmo, em uma máquina com $M$ linhas de cache ideal, cada uma com tamanho $b$ bytes, é limitado por

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

Algoritmos sub-cúbicos

Melhoria das estimativas do expoente

ω ao

longo do tempo para a complexidade computacional da multiplicação de matrizes .

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

Existem algoritmos que fornecem melhores tempos de execução do que os simples. O primeiro a ser descoberto foi o algoritmo de Strassen , desenvolvido por Volker Strassen em 1969 e frequentemente referido como "multiplicação rápida da matriz". É baseado em uma maneira de multiplicar duas matrizes $2 \times 2$ que requer apenas 7 multiplicações (em vez das 8 usuais), à custa de várias operações adicionais de adição e subtração. Aplicar isso recursivamente fornece um algoritmo com um custo multiplicativo de . O algoritmo de Strassen é mais complexo e a estabilidade numérica é reduzida em comparação ao algoritmo ingênuo, mas é mais rápido nos casos em que $n$ $> 100$ ou mais e aparece em várias bibliotecas, como BLAS . É muito útil para grandes matrizes sobre domínios exatos, como campos finitos , onde a estabilidade numérica não é um problema. ${\ displaystyle O (n ^ {\ log _ {2} 7}) \ aprox O (n ^ {2.807})}$

É uma questão em aberto na ciência da computação teórica quão bem o algoritmo de Strassen pode ser melhorado em termos de complexidade assintótica . O expoente de multiplicação da matriz , geralmente denotado , é o menor número real para o qual qualquer matriz sobre um campo pode ser multiplicada em conjunto usando operações de campo. O melhor limite atual é , por Josh Alman e Virginia Vassilevska Williams . Este algoritmo, como todos os outros algoritmos recentes nesta linha de pesquisa, é uma generalização do algoritmo Coppersmith – Winograd, que foi dado por Don Coppersmith e Shmuel Winograd em 1990. A ideia conceitual desses algoritmos é semelhante ao algoritmo de Strassen: um caminho é desenvolvido para multiplicar duas matrizes $k$ $\times$ $k$ com menos de $k$ $3$ multiplicações, e esta técnica é aplicada recursivamente. No entanto, o coeficiente constante oculto pela notação Big O é tão grande que esses algoritmos só valem a pena para matrizes muito grandes para serem manipuladas nos computadores atuais. ${\ displaystyle \ omega}$ ${\ displaystyle n \ times n}$ ${\ displaystyle n ^ {\ omega + o (1)}}$ ${\ displaystyle \ omega}$ ${\ displaystyle \ omega <2,3728596}$

Algoritmo Freivalds é um simples algoritmo de Monte Carlo que, matrizes dadas $Um$ , $B$ e $C$ , verifica em $Θ (n 2)$ tempo se $AB = C$ .

Algoritmos paralelos e distribuídos

Paralelismo de memória compartilhada

O algoritmo de divisão e conquista esboçado anteriormente pode ser paralelizado de duas maneiras para multiprocessadores de memória compartilhada . Estes são baseados no fato de que as oito multiplicações recursivas de matrizes em

{\ 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} e A_ {21} B_ {12} + A_ {22} B_ {22} \\\ end {pmatriz}}}

podem ser executados independentemente um do outro, assim como as quatro somas (embora o algoritmo precise "juntar" as multiplicações antes de fazer as somas). Explorando todo o paralelismo do problema, obtém-se um algoritmo que pode ser expresso em pseudocódigo de estilo fork – join :

$Multiplicação do$ procedimento $($ $C$ $,$ $A$ $,$ $B$ $)$ :

Caso base: se $n = 1$ , defina $c 11 \leftarrow a 11 \times b 11$ (ou multiplique uma matriz de bloco pequeno).
Caso contrário, aloque espaço para uma nova matriz T de forma n × n , então:
- Divida $A$ em $A 11$ , $A 12$ , $A 21$ , $A 22$ .
- Divida $B$ em $B 11$ , $B 12$ , $B 21$ , $B 22$ .
- Particionar $C$ em $C 11$ , $C 12$ , $C 21$ , $C 22$ .
- Divida $T$ em $T 11$ , $T 12$ , $T 21$ , $T 22$ .
- Execução paralela:
  - $Multiplique o$ garfo $($ $C$ $11$ $,$ $A$ $11$ $,$ $B$ $11$ $)$ .
  - $Multiplique o$ garfo $($ $C$ $12$ $,$ $A$ $11$ $,$ $B$ $12$ $)$ .
  - $Multiplique o$ garfo $($ $C$ $21$ $,$ $A$ $21$ $,$ $B$ $11$ $)$ .
  - $Multiplique o$ garfo $($ $C$ $22$ $,$ $A$ $21$ $,$ $B$ $12$ $)$ .
  - $Multiplique o$ garfo $($ $T$ $11$ $,$ $A$ $12$ $,$ $B$ $21$ $)$ .
  - $Multiplique o$ garfo $($ $T$ $12$ $,$ $A$ $12$ $,$ $B$ $22$ $)$ .
  - $Multiplique o$ garfo $($ $T$ $21$ $,$ $A$ $22$ $,$ $B$ $21$ $)$ .
  - $Multiplique o$ garfo $($ $T$ $22$ $,$ $A$ $22$ $,$ $B$ $22$ $)$ .
- Junte-se (aguarde a conclusão dos garfos paralelos).
- $adicione (C, T)$ .
- Desalocar $T$ .

O procedimento $add (C, T)$ adiciona $T$ em $C$ , elemento a elemento:

Caso base: se $n = 1$ , defina $c 11 \leftarrow c 11 + t 11$ (ou faça um loop curto, talvez desenrolado).
De outra forma:
- Particionar $C$ em $C 11$ , $C 12$ , $C 21$ , $C 22$ .
- Divida $T$ em $T 11$ , $T 12$ , $T 21$ , $T 22$ .
- Em paralelo:
  - $Adicionar$ garfo $($ $C$ $11$ $,$ $T$ $11$ $)$ .
  - $Adicionar$ garfo $($ $C$ $12$ $,$ $T$ $12$ $)$ .
  - $Adicionar$ garfo $($ $C$ $21$ $,$ $T$ $21$ $)$ .
  - $Adicionar$ garfo $($ $C$ $22$ $,$ $T$ $22$ $)$ .
- Junte-se .

Aqui, fork é uma palavra-chave que sinaliza que um cálculo pode ser executado em paralelo com o resto da chamada de função, enquanto join aguarda a conclusão de todos os cálculos previamente "bifurcados". $partição$ atinge seu objetivo apenas pela manipulação do ponteiro.

Este algoritmo tem um comprimento de caminho crítico de $Θ (log 2 n)$ etapas, o que significa que leva muito tempo em uma máquina ideal com um número infinito de processadores; portanto, tem uma aceleração máxima possível de $Θ (n 3 / log 2 n)$ em qualquer computador real. O algoritmo não é prático devido ao custo de comunicação inerente à movimentação de dados de e para a matriz temporária $T$ , mas uma variante mais prática atinge $Θ (n 2)$ speedup, sem usar uma matriz temporária.

Multiplicação de matrizes de blocos. No algoritmo 2D, cada processador é responsável por uma submatriz de

C

. No algoritmo 3D, cada par de submatrizes de

A

e

B

que é multiplicado é atribuído a um processador.

Algoritmos que evitam a comunicação e distribuídos

Em arquiteturas modernas com memória hierárquica, o custo de carregamento e armazenamento de elementos da matriz de entrada tende a dominar o custo da aritmética. Em uma única máquina, essa é a quantidade de dados transferidos entre a RAM e o cache, enquanto em uma máquina com vários nós de memória distribuída é a quantidade transferida entre os nós; em ambos os casos, é chamada de largura de banda de comunicação . O algoritmo ingênuo usando três loops aninhados usa largura de banda de comunicação $Ω (n 3)$ .

O algoritmo de Cannon , também conhecido como algoritmo 2D , é um algoritmo que evita a comunicação que particiona cada matriz de entrada em uma matriz de bloco cujos elementos são submatrizes de tamanho $\sqrt M / 3$ por $\sqrt M / 3$ , onde $M$ é o tamanho da memória rápida. O algoritmo ingênuo é então usado sobre as matrizes de bloco, computando produtos de submatrizes inteiramente na memória rápida. Isso reduz a largura de banda de comunicação para $O (n 3 / \sqrt M)$ , que é assintoticamente ideal (para algoritmos realizando computação $Ω (n 3)$ ).

Em uma configuração distribuída com $p$ processadores dispostos em uma malha $\sqrt p$ por $\sqrt p$ 2D, uma submatriz do resultado pode ser atribuída a cada processador e o produto pode ser calculado com cada processador transmitindo $O (n 2 / \sqrt p)$ palavras, que é assintoticamente ideal assumindo que cada nó armazena o mínimo de $O (n 2 / p)$ elementos. Isso pode ser melhorado pelo algoritmo 3D, que organiza os processadores em uma malha de cubo 3D, atribuindo cada produto de duas submatrizes de entrada a um único processador. As submatrizes de resultado são geradas realizando uma redução em cada linha. Este algoritmo transmite $O (n 2 / p 2/3)$ palavras por processador, que é assintoticamente ideal. No entanto, isso requer a replicação de cada elemento da matriz de entrada $p 1/3$ vezes e, portanto, requer um fator de $p 1/3 a$ mais de memória do que o necessário para armazenar as entradas. Este algoritmo pode ser combinado com Strassen para reduzir ainda mais o tempo de execução. Os algoritmos "2.5D" fornecem uma compensação contínua entre o uso de memória e a largura de banda de comunicação. Em ambientes modernos de computação distribuída, como MapReduce , algoritmos de multiplicação especializados foram desenvolvidos.

Algoritmos para malhas

Multiplicação de matriz concluída em etapas 2n-1 para duas matrizes n × n em uma malha de arame cruzado.

Existem vários algoritmos para multiplicação em malhas . Para a multiplicação de dois n × n em uma malha bidimensional padrão usando o algoritmo do Canhão 2D , pode-se completar a multiplicação em 3 n -2 etapas, embora isso seja reduzido à metade desse número para cálculos repetidos. A matriz padrão é ineficiente porque os dados das duas matrizes não chegam simultaneamente e devem ser preenchidos com zeros.

O resultado é ainda mais rápido em uma malha com fiação cruzada de duas camadas, onde apenas 2 n -1 etapas são necessárias. O desempenho melhora ainda mais para cálculos repetidos, levando a uma eficiência de 100%. A matriz de malha com fio cruzado pode ser vista como um caso especial de uma estrutura de processamento não plana (ou seja, multicamadas).

Veja também

Referências

Leitura adicional

Buttari, Alfredo; Langou, Julien; Kurzak, Jakub; Dongarra, Jack (2009). "Uma classe de algoritmos paralelos de álgebra linear para arquiteturas multicore". Computação 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). "Anatomia da multiplicação de matrizes de alto desempenho". Transações ACM em 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: uma estrutura para instanciar rapidamente a funcionalidade do BLAS". Transações ACM em software matemático . 41 (3): 1–33. doi : 10.1145 / 2764454 . S2CID 1242360 .
Como otimizar o GEMM

Languages

In other projects

Algoritmo de multiplicação de matriz - Matrix multiplication algorithm

Conteúdo

Algoritmo iterativo

Comportamento do cache

Algoritmo de divisão e conquista

Matrizes não quadradas

Comportamento do cache

Algoritmos sub-cúbicos

Algoritmos paralelos e distribuídos

Paralelismo de memória compartilhada

Algoritmos que evitam a comunicação e distribuídos

Algoritmos para malhas

Veja também

Referências

Leitura adicional