O produto de Kronecker é um tipo especial de produto de matrizes (produto tensorial).

Esse produto é aplicado, por exemplo, no cálculo do laplaciano discreto, que, por sua vez, é utilizado na resolução numérica de algumas equações diferenciais parciais, como a equação de Laplace.
O objetivo desta postagem é apresentar um algoritmo que permite calcular o produto de Kronecker entre duas matrizes.
Além do pseudocódigo, disponibilizarei o algoritmo codificado em Java e em C.
Também mostrarei como calcular o produto de Kronecker no MATLAB e no Octave e um exemplo de equação resolvida numericamente com o auxílio desse produto.
Produto de Kronecker
Se Ap×q é uma matriz com p linhas e q colunas e Br×s é uma matriz com r linhas e s colunas, então o produto de Kronecker de A e B (denotado por A⊗B) é uma matriz Cpr×qs com pr linhas e qs colunas dividida em p×q blocos. Cada bloco (i,j) é uma matriz da forma aijB, onde aij é o elemento da matriz A que ocupa a linha i e coluna j (TREFETHEN, 2000).
Não entendeu nada? Um exemplo deixará mais claro:

Em geral, o produto de Kronecker não é comutativo
A⊗B≠B⊗A
Além disso, não há restrições em relação às dimensões de A e de B.
[24]⊗[2131]=[426284124]
A maneira como os blocos são distribuídos depende da matriz A. Nesse último exemplo, a matriz A tem tamanho 2×1. Portanto, o produto conterá 2×1 blocos. O primeiro bloco (1, 1) é igual ao elemento 2 multiplicado pela matriz B. O segundo bloco (2, 1) é igual ao produto do elemento 4 por B.
Na prática, o produto de Kronecker é mais fácil do que o produto usual de matrizes. Por outro lado, seu algoritmo é um pouco mais complexo.
Os blocos mencionados ao longo do texto são apenas uma maneira de facilitar a compreensão, pois o resultado do produto A⊗B é uma matriz como qualquer outra.
Algoritmo
O algoritmo do produto de Kronecker é apresentado a seguir
1. kron(Ap×q, Br×s) 2. | inicializar a matriz Cpr×qs 3. | para i ← 1 até p 4. | | para j ← 1 até q 5. | | | row ← (i - 1) * r + 1 6. | | | para k ← até r 7. | | | | col ← (j - 1) * s + 1 8. | | | | para l ← 1 até s 9. | | | | | C(row, col) ← A(i, j) * B(k, l) 10.| | | | | col ← col + 1 11.| | | | fim_para 12.| | | | row ← row + 1 13.| | | fim_para 14.| | fim_para 15.| fim_para 16.| retorne Cpr×qs 17.fim_kron
A complexidade no tempo desse algoritmo é O(p×r×q×s).
Códigos
Os códigos do produto de Kronecker em Java e em C são apresentados a seguir.
Código em Java
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 | //GitHub: HenriqueIni / www.blogcyberini.com public class ProdutoKronecker { //Calcula o produto de Kronecker entre A e B public static double [][] kron( double [][] A, double [][] B){ //O código não faz validações. A e B devem ser matrizes válidas int p = A.length; int q = A[ 0 ].length; int r = B.length; int s = B[ 0 ].length; int row, col; //Matriz resultado double [][] C = new double [p * r][q * s]; //Laço quádruplo aninhado: percorre as matrizes A e B for ( int i = 0 ; i < p; i++){ for ( int j = 0 ; j < q; j++){ row = i * r; for ( int k = 0 ; k < r; k++){ col = j * s; for ( int l = 0 ; l < s; l++){ //Preenche a matriz C C[row][col] = A[i][j] * B[k][l]; col++; } row++; } } } return C; } // Método auxiliar para transformar matriz em String private static String toString( double [][] matrix){ String aux = "" ; for ( int i = 0 ; i < matrix.length; i++){ aux += "[" + matrix[i][ 0 ]; for ( int j = 1 ; j < matrix[i].length; j++){ aux += "," + matrix[i][j]; } aux += "]" ; if (i != matrix.length - 1 ){ aux += "\n" ; } } return aux; } //Testes public static void main(String[] args) { //Teste 1 double [][] A = {{ 2 , 4 }}; double [][] B = {{ 2 , 1 },{ 3 , 1 }}; System.out.println(toString(kron(A, B))+ "\n" ); //Teste 2 double [][] C = {{ 1 , 2 }}; double [][] D = {{ 1 , 3 }}; double [][] E = {{ 1 },{ 1 }}; System.out.println(toString(kron(kron(C, D), E))+ "\n" ); //Teste 3 double [][] F = {{ 2 },{ 4 }}; System.out.println(toString(kron(F, B))+ "\n" ); } } |
Código em C
Observe que no código em C a matriz que vai conter os resultados deve ser passada nos parâmetros. Se você for mais familiarizado com ponteiros em C (não é o meu caso), então você pode deixar esse algoritmo mais simples ainda.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | //GitHub: HenriqueIni / www.blogcyberini.com #include <stdio.h> #include <stdlib.h> //Calcula o produto de Kronecker entre A e B //O resultado é armazenado em C //Esse método pode ser simplificado com o uso de ponteiros void kron( int p, int q, double A[p][q], int r, int s, double B[r][s], int pr, int qs, double C[pr][qs]){ //O código não faz validações, logo os parâmetros devem ser válidos int i, j, k, l, col, row; for (i = 0; i < p; i++){ for (j = 0; j < q; j++){ row = i * r; for (k = 0; k < r; k++){ col = j * s; for (l = 0; l < s; l++){ //Preenche a matriz C C[row][col] = A[i][j] * B[k][l]; col++; } row++; } } } } // Função auxiliar para imprimir matrizes void printArray( int n, int m, double A[n][m]){ int i, j; for (i = 0; i < n; i++){ printf ( "[%.4f" , A[i][0]); for (j = 1; j < m; j++){ printf ( ",%.4f" , A[i][j]); } printf ( "]" ); if (i != n - 1){ printf ( "\n" ); } } } // Testes int main() { //Teste 1 double A[1][2] = {{2,4}}; double B[2][2] = {{2, 1},{3, 1}}; double R[2][4]; kron(1, 2, A, 2, 2, B, 2, 4, R); printArray(2, 4, R); printf ( "\n\n" ); //Teste 2 double F[2][1] = {{2},{4}}; double R2[4][2]; kron(2, 1, F, 2, 2, B, 4, 2, R2); printArray(4, 2, R2); } |
MATLAB e Octave
Os softwares MATLAB e Octave já possuem uma função para calcular o produto de Kronecker.

Em ambos os programas, o produto de Kronecker entre as matrizes A
e B
pode ser feito utilizando o seguinte trecho de código [2][3]:
C = kron(A, B)
Laplaciano discreto
Se D2n é uma matriz de derivação n×n de segunda ordem em uma dimensão (isto é, ela serve para calcular a derivada de segunda ordem de uma função de uma variável), então o laplaciano discreto em duas dimensões será [1]
Ln=In⊗D2n+D2n⊗In
In é a matriz identidade de ordem n. Em três dimensões, teríamos
Ln=In⊗In⊗D2n+In⊗D2n⊗In+D2n⊗In⊗In
Nos dois exemplos, estou supondo que uma mesma matriz de derivação é utilizada para derivar em todas as dimensões.
Aplicação
O exemplo a seguir foi retirado do livro "Spectral Methods in MATLAB" por Trefethen (2000).
No exemplo, Trefethen resolve numericamente a equação de Poisson ∇2u=10sin(8x(y−1)) no retângulo [−1,1]×[−1,1], onde a função u é igual a zero na fronteira.
Para resolver a equação, o laplaciano discreto é calculado utilizando o produto de Kronecker.

O programa utilizado nos cálculos pode ser obtido no link: Equação de Poisson (Octave). Eu fiz algumas modificações para que o programa funcionasse no Octave, uma vez que o original foi escrito para MATLAB.
Se quiser uma visão mais ampla desse gráfico, assista ao vídeo: gráfico da solução da equação (YouTube).
Leia também
Referências
- [1] Trefethen, Lloyd. Spectral Methods in MATLAB. SIAM, 2000.
- [2] Octave DOCS: 18.4 Functions of a Matrix. Acesso em 16 de abril de 2018.
- [3] MathWorks: kron. Acesso em 16 de abril de 2018.
Nenhum comentário:
Postar um comentário