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 $A_{p\times q}$ é uma matriz com $p$ linhas e $q$ colunas e $B_{r\times s}$ é uma matriz com $r$ linhas e $s$ colunas, então o produto de Kronecker de $A$ e $B$ (denotado por $A\otimes B$) é uma matriz $C_{pr\times qs}$ com $pr$ linhas e $qs$ colunas dividida em $p\times q$ blocos. Cada bloco $(i, j)$ é uma matriz da forma $a_{ij} B$, onde $a_{ij}$ é 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

$$A\otimes B\neq B\otimes A$$

Além disso, não há restrições em relação às dimensões de A e de B.

$$\begin{bmatrix}2\\4\end{bmatrix}\otimes\begin{bmatrix}2 & 1\\3 & 1\end{bmatrix}=\begin{bmatrix}4 & 2\\6 & 2\\\hline8 & 4\\12 & 4\end{bmatrix}$$

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\otimes 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\times r\times q\times s)$.

Códigos

Os códigos do produto de Kronecker em Java e em C são apresentados a seguir.

Código em Java

//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.

//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 $D_n^2$ é uma matriz de derivação $n\times 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]

$$L_n=I_n\otimes D_n^2 + D_n^2\otimes I_n$$

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

$$L_n=I_n\otimes I_n\otimes D_n^2 + I_n\otimes D_n^2\otimes I_n + D_n^2\otimes I_n\otimes I_n$$

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 $\nabla^2 u = 10\operatorname{sin}(8x(y - 1))$ no retângulo $[-1,1]\times[-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


Siga o blog

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

Receba notificações de novas postagens e novidades do blog por e-mail.

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