Algoritmo do Produto de Kronecker

Por
|  

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

Algoritmo do Produto de Kronecker

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 AB) é 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:

Exemplo de produto de Kronecker
Exemplo de produto de Kronecker [1].

Em geral, o produto de Kronecker não é comutativo

ABBA

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 AB é 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.

Produto de Kronecker no Octave
Exemplo de produto de Kronecker no Octave

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=InD2n+D2nIn

In é a matriz identidade de ordem n. Em três dimensões, teríamos

Ln=InInD2n+InD2nIn+D2nInIn

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(y1)) 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.

Solução da equação de Poisson no Octave
Exemplo de solução da equação de Poisson no Octave [1]

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

Sugestões de livros para estudantes de computação na Amazon (patrocinado): Lista de Livros

Obrigado pela leitura! Se você puder, considere apoiar financeiramente o Blog Cyberini, Chave Pix: cyberpix9@gmail.com

Doar com PayPal

Siga o blog

Redes sociais: Facebook, Twitter, YouTube, Pinterest, Instagram, Telegram

Importante: utilize o bom senso na hora de comentar. Acesse a política de privacidade para maiores informações sobre comentários.

Nenhum comentário:

Postar um comentário

Este site usa cookies para fornecer serviços, analisar o tráfego e exibir publicidade personalizada, conforme a nossa política de privacidade. Seus dados, como IP e user agent, são compartilhados com nossos parceiros (como o Google). Saiba Mais

RecusarAceitar