Neste artigo, apresento um algoritmo que permite calcular determinantes em tempo O(n3), sendo mais eficiente que o teorema de Laplace, cuja complexidade é O(n!).
Índice
Triangularização de uma Matriz
O algoritmo proposto consiste em obter uma matriz triangular com o mesmo determinante da matriz na qual desejamos calcular o determinante.
Isso é possível graças às seguintes propriedades
Propriedade 1. Troca de filas paralelas: se trocarmos de posição duas linhas (ou duas colunas) de uma matriz quadrada M, o determinante da nova matriz obtida é o oposto do determinante da matriz anterior (DANTE, 2015).
Propriedade 2. Teorema de Jacobi: Seja A uma matriz quadrada. Se multiplicarmos todos os elementos de uma linha (ou coluna) pelo mesmo número e somarmos os resultados aos elementos correspondentes de outra linha (ou coluna), formando a matriz B, então det A = det B (DANTE, 2015).
Para exemplificá-las, vamos utilizar a matriz a seguir
$$A=\begin{bmatrix}1 & 5\\-2 & 3\end{bmatrix}.$$
O determinante de A é igual a 13. Se trocarmos a primeira coluna pela segunda, obteremos a seguinte matriz
$$B=\begin{bmatrix}5 & 1\\3 & -2\end{bmatrix},$$
cujo determinante é igual a -13. Ou seja, o oposto de det A.
Para exemplificar a segunda propriedade, vamos multiplicar a primeira linha de A por 2 e somar os resultados à segunda linha e, dessa forma, obteremos a seguinte matriz
$$C=\begin{bmatrix}1 & 5\\0 & 13\end{bmatrix}.$$
Observe que det C = 13, isto é, det C = det A. Além disso, a nossa escolha foi bem conveniente, pois C é uma matriz triangular.
O determinante de uma matriz triangular é igual ao produto dos elementos da diagonal principal e isso pode ser realizado em tempo O(n).
O algoritmo de triangularização é análogo à eliminação de Gauss, utilizada para escalonar sistemas lineares. Optaremos por um algoritmo que calcula uma matriz triangular superior B a partir de uma matriz quadrada A tal que $\det B = \det A$. Tenha em mente que é uma questão de escolha e que B poderia ser uma matriz triangular inferior, entretanto o algoritmo seria diferente.
A ideia do algoritmo é utilizar o teorema de Jacobi para "criar os zeros" necessários para que a matriz A se torne a matriz triangular B. Um exemplo deixará isso mais claro. Seja a matriz 3×3 a seguir
$$\begin{bmatrix}2 & -2 & -1\\3 & -4 & 1\\1 & 1 & 5\end{bmatrix}.$$
Vamos "zerar" (não achei termo melhor) os elementos da primeira coluna nas linhas 2 e 3 utilizando a linha 1. No caso da linha 2, temos que multiplicar a linha 1 por um fator F tal que 2F = -3, onde 2 é o elemento a11 e -3 é o oposto de a21. Neste caso, F = -3/2, logo
$$\begin{bmatrix}2 & -2 & -1\\3 + 2F & -4 + (-2)\times F & 1+(-1)\times F\\1 & 1 & 5\end{bmatrix}.$$
$$\begin{bmatrix}2 & -2 & -1\\3 + (-3) & -4 + 3 & 1+3/2\\1 & 1 & 5\end{bmatrix}.$$
$$\begin{bmatrix}2 & -2 & -1\\0 & -1 & 5/2\\1 & 1 & 5\end{bmatrix}.$$
Para a linha 3, o fator será F = -1/2
$$\begin{bmatrix}2 & -2 & -1\\0 & -1 & 5/2\\1 + (-1) & 1+1 & 5+1/2\end{bmatrix}.$$
$$\begin{bmatrix}2 & -2 & -1\\0 & -1 & 5/2\\0 & 2 & 11/2\end{bmatrix}.$$
Agora que já utilizamos a linha 1, vamos utilizar a linha 2 para "zerar" todos os elementos da coluna 2 nas linhas posteriores à linha 2 (neste caso, apenas a linha 3). O fator F será F = 2
$$\begin{bmatrix}2 & -2 & -1\\0 & -1 & 5/2\\0 & 2+(-2) & 11/2 + 5\end{bmatrix}.$$
$$\begin{bmatrix}2 & -2 & -1\\0 & -1 & 5/2\\0 & 0 & 21/2\end{bmatrix}.$$
Isso conclui a triangularização. Obtemos o determinante multiplicando os elementos da diagonal principal, que é igual a -21.
Algoritmo de Triangularização
A triangularização que exemplificamos é expressa pelo seguinte algoritmo:
1. triang(A, n) //n é a ordem da matriz 2. para k ← 1 até n - 1 3. para m ← k + 1 até n 4. F ← (-1) * A(m, k)/A(k, k) 5. para l ← k até n 6. A(m, l) ← A(m, l) + F * A(k, l) 7. fim_para 8. fim_para 9. fim_para 10.fim_triang
Observe que o algoritmo triang(A, n) apenas realiza a triangularização e que a indexação inicia por 1. A notação A(i, j) indica o elemento da linha i e coluna j.
O loop na linha 2 percorre cada linha (exceto a última). A k-ésima linha é utilizada para "criar zeros" na coluna k, cujas linhas estão entre k+1 e n (loop da linha 3). A linha 4 calcula o fator F, que é o mesmo do exemplo dado anteriormente. O loop das linhas 5-7 multiplica os elementos da linha k por F e soma os resultados aos elementos correspondentes da linha m.
Observe que só precisamos considerar os elementos da coluna k em diante, pois os elementos anteriores já são iguais a zero (eles foram "zerados" em iterações anteriores). Observe também que a primeira iteração em 5-7 é irrelevante, pois ela basicamente faz A(m, k) ← 0, ou seja, você pode eliminar essa iteração e fazer a atribuição diretamente.
O algoritmo apresenta um erro: se A(k, k) for igual a zero (linha 4) teríamos uma divisão por zero! Tal problema pode ser resolvido aplicando a primeira propriedade apresentada, isto é, trocamos a k-ésima linha por outra linha cujo elemento da coluna k seja diferente de zero. É claro, não pode ser qualquer linha, pois devemos manter os zeros criados nas iterações anteriores em suas devidas posições. Para que isso ocorra, devemos escolher uma linha posterior à linha k.
Se um elemento diferente de zero não for encontrado na coluna k, a partir da linha k+1, então o determinante será zero. Por que zero? Simples! É porque a submatriz dos elementos das linhas i≥k e colunas j≥k tem determinante igual a zero, pois a primeira coluna da submatriz só tem elementos iguais a zero. Um exemplo ilustrará a situação. Suponha que, durante a execução do algoritmo numa matriz 5×5, a matriz assuma a seguinte forma
A área em destaque é a submatriz com uma fila de zeros. Se aplicarmos o teorema de Laplace para o cálculo do determinante, obteremos
$$a_{11}\times a_{22}\times\begin{bmatrix}0 & a_{34} & a_{35}\\0 & a_{44} & a_{45}\\0 & a_{54} & a_{55}\end{bmatrix}$$
Que obviamente é igual a zero, pois a submatriz tem uma fila de zeros. Aplicando as correções, obtemos o seguinte algoritmo
1. triang(A, n) 2. | p ← 1 3. | para k ← 1 até n - 1 4. | | se(A(k, k) = 0) 5. | | | para j ← k + 1 até n 6. | | | | se(A(j, k) ≠ 0) 7. | | | | | trocar linhas j e k 8. | | | | | p ← p * (-1) 9. | | | | | encerrar loop (linha 5) 10. | | | | fim_se 11. | | | fim_para 12. | | fim_se 13. | | se(A(k, k) = 0) 14. | | | A é uma matriz singular (det A = 0) 15. | | fim_se 16. | | para m ← k + 1 até n 17. | | | F ← (-1) * A(m, k)/A(k, k) 18. | | | para l ← k até n 19. | | | | A(m, l) ← A(m, l) + F * A(k, l) 20. | | | fim_para 21. | | fim_para 22. | fim_para 23. fim_triang
A variável p no pseudocódigo serve para ajustar o sinal do determinante. O algoritmo só realiza a triangularização. Após calcular o determinante da matriz triangular, é necessário multiplicar o determinante por p. Apesar de estar matematicamente correto, esse algoritmo ainda possui um problema que discutiremos na seção a seguir.
Erros de Arredondamento
O problema do algoritmo apresentado anteriormente está na divisão A(m, k)/A(k, k), principalmente quando A(k, k) é um valor próximo de zero. A representação em ponto flutuante nos computadores não permite a representação exata de números decimais infinitos e isso gera erros de arredondamento que se propagam ao longo dos cálculos.
Para amenizar esses erros, podemos utilizar a propriedade da troca de linhas de tal forma que A(k, k) seja sempre o maior valor possível em módulo. Isso ajuda a reduzir os erros, mas não resolve completamente o problema. Dependendo do quão mal condicionada a matriz for, os resultados podem ser insatisfatórios e outro algoritmo deve ser utilizado. Para maiores informações sobre o assunto, consulte um livro de Cálculo Numérico (deixarei um nas referências). O algoritmo final é apresentado a seguir
1. detTriang(A, n) 2. | p ← 1 //fator de ajuste do sinal 3. | para k ← 1 até n - 1 4. | | max ← |A(k ,k)| 5. | | maxIndex ← k 6. | | para i ← k + 1 até n 7. | | | se(max < |A(i, k)|) 8. | | | | max ← |A(i, k)| 9. | | | | maxIndex ← i 10. | | | fim_se 11. | | fim_para 12. | | se(maxIndex ≠ k) 13. | | | p ← p * (-1) 14. | | | para j ← 1 até n 15. | | | | temp ← A(k, j) 16. | | | | A(k, j) ← A(maxIndex, j) 17. | | | | A(maxIndex, j) ← temp 18. | | | fim_para 19. | | fim_se 20. | | se(A(k, k) = 0) 21. | | | retorne 0 //a matriz é singular 22. | | senão 23. | | | para m ← k + 1 até n 24. | | | | F ← (-1) * A(m, k)/A(k, k) 25. | | | | A(m, k) ← 0 //elimina a primeira iteração 26. | | | | para l ← k + 1 até n 27. | | | | | A(m, l) ← A(m, l) + F * A(k, l) 28. | | | | fim_para 29. | | | fim_para 30. | | fim_se 31. | fim_para 32. | //calcula o determinante da matriz 33. | det ← 1 34. | para q ← 1 até n 35. | | det ← det * A(q, q) 36. | fim_para 37. | retorne p * det 38. fim_triang
As linhas 4-19 encontram o maior valor em módulo na coluna k, a partir da linha k, fazendo a troca da linha k pela linha maxIndex, que contém o valor encontrado. Se a linha com o maior valor em módulo for k, então a troca não é realizada, obviamente.
Análise do Algoritmo
O algoritmo apresentado tem complexidade no tempo de O(n3), conforme o triplo loop aninhado sugere. Considerando que o teorema de Laplace tem complexidade O(n!), o algoritmo é muito mais eficiente. Além disso, não é necessário alocar memória extra e os resultados dos cálculos podem ser armazenados na própria matriz (ou numa cópia, caso você precise manter a matriz original). Já o teorema de Laplace requer alocação de memória extra para gerar as matrizes do cálculo dos cofatores, algo que representa um grande custo no espaço.
Códigos
A seguir, apresento o algoritmo codificado nas linguagens Java e C. Os códigos modificam a matriz passada como parâmetro, convertendo-a na matriz triangular utilizada no cálculo do determinante. Portanto, se for necessário manter a matriz original, você deve passar como parâmetro uma cópia da matriz ou modificar o algoritmo para que ele crie uma cópia da matriz e realize os cálculos com a cópia (o algoritmo continuaria com complexidade O(n3)). Além disso, os algoritmos não realizam validações, logo os parâmetros passados para os algoritmos devem ser válidos.
Código em Java
Ao invés de passar a ordem da matriz como parâmetro, você pode obter a ordem da seguinte forma: int n = A.length. Obviamente, isso não exclui a necessidade da matriz precisar ser válida.
//Código por Henrique Felipe (GitHub: HenriqueIni) public class DetTriang { //Calcula o determinante da matriz A, de ordem n, via triangularização [O(n³)]. public static double detTriang(double[][] A, int n) { //O código não faz validações, portanto a matriz precisa ser válida int p = 1; //fator de ajuste do sinal for (int k = 0; k < n - 1; k++) { double max = Math.abs(A[k][k]); int maxIndex = k; //encontra o maior valor em módulo na coluna k, a partir da linha k for (int i = k + 1; i < n; i++) { if (max < Math.abs(A[i][k])) { max = Math.abs(A[i][k]); maxIndex = i; } } //troca a linha k pela linha com o maior valor absoluto na coluna k if (maxIndex != k) { double temp; p = p * (-1); for (int j = 0; j < n; j++) { temp = A[k][j]; A[k][j] = A[maxIndex][j]; A[maxIndex][j] = temp; } } //se depois de tudo A[k][k] for zero, então a matriz é singular if (A[k][k] == 0) { return 0; } else { //'cria zeros' na coluna k, depois da linha k for (int m = k + 1; m < n; m++) { double F = (-1) * A[m][k] / A[k][k]; A[m][k] = 0; //evita a primeira iteração for (int l = k + 1; l < n; l++) { A[m][l] = A[m][l] + F * A[k][l]; } } } } //depois da triangularização, calcula o determinante double det = 1; for (int q = 0; q < n; q++) { //multiplica os elementos na diagonal principal det = det * A[q][q]; } return p * det; } //Testes public static void main(String[] args) { double[][] A = {{2, -2, -1}, {3, -4, 1}, {1, 1, 5}}; //det A = -21 double[][] B = {{1, 5}, {-2, 3}}; //det B = 13 double[][] C = {{1, 2, 3, 4}, {4, 3, 2, 1}, {5, 6, 7, 8}, {8, 7, 6, 5}}; //det C = 0 double[][] D = {{-2, 3, 1, 7}, {0, -1, 2, 1}, {3, -4, 5, 1}, {1, 0, -2, -1}}; //det D = 24 double[][] E = {{1, 0, 5, 0}, {2, -1, 0, 3}, {3, 0, 2, 0}, {7, 0, 6, 5}}; //det E = 65 double[][] F = {{1, 2, 3, 4, 5}, {-1, 0, 8, 5, 9}, {2, 2, 3, 5, 7}, {-3, -1, 12, 0, 7}, {1, -1, 1, 3, 13}}; //det F = -710 System.out.println("det A = " + detTriang(A, 3)); System.out.println("det B = " + detTriang(B, 2)); System.out.println("det C = " + detTriang(C, 4)); System.out.println("det D = " + detTriang(D, 4)); System.out.println("det E = " + detTriang(E, 4)); System.out.println("det F = " + detTriang(F, 5)); } }
Código em C
//Código por Henrique Felipe (GitHub: HenriqueIni) #include <stdio.h> #include <stdlib.h> #include <math.h> //Calcula o determinante da matriz A, de ordem n, via triangularização [O(n³)]. double detTriang(int n, double A[n][n]){ //O código não faz validações, portanto a matriz precisa ser válida int p = 1; //fator de ajuste do sinal int k; for(k = 0; k < n - 1; k++){ double max = abs(A[k][k]); int maxIndex = k; //encontra o maior valor em módulo na coluna k, a partir da linha k int i; for(i = k + 1; i < n; i++){ if(max < abs(A[i][k])){ max = abs(A[i][k]); maxIndex = i; } } //troca a linha k pela linha com o maior valor absoluto na coluna k if(maxIndex != k){ double temp; p = p * (-1); int j; for(j = 0; j < n; j++){ temp = A[k][j]; A[k][j] = A[maxIndex][j]; A[maxIndex][j] = temp; } } //se depois de tudo A[k][k] for zero, então a matriz é singular if(A[k][k] == 0){ return 0; }else{ //'cria zeros' na coluna k, depois da linha k int m; for(m = k + 1; m < n; m++){ double F = (-1) * A[m][k]/A[k][k]; A[m][k] = 0; //evita a primeira iteração int l; for(l = k + 1; l < n; l ++){ A[m][l] = A[m][l] + F * A[k][l]; } } } } //depois da triangularização, calcula o determinante double det = 1; int q; for(q = 0; q < n; q++){ //multiplica os elementos na diagonal principal det = det * A[q][q]; } return p * det; } //Testes int main() { double A[3][3] = { {2, -2, -1}, {3, -4, 1}, {1, 1, 5}}; //det A = -21 double B[2][2] = { {1, 5}, {-2, 3}}; //det B = 13 double C[4][4] = { {1, 2, 3, 4}, {4, 3, 2, 1}, {5, 6, 7, 8}, {8, 7, 6, 5}}; //det C = 0 double D[4][4] = { {-2, 3, 1, 7}, {0, -1, 2, 1}, {3, -4, 5, 1}, {1, 0, -2, -1}}; //det D = 24 double E[4][4] = { {1, 0, 5, 0}, {2, -1, 0, 3}, {3, 0, 2, 0}, {7, 0, 6, 5}}; //det E = 65 double F[5][5] = { {1, 2, 3, 4, 5}, {-1, 0, 8, 5, 9}, {2, 2, 3, 5, 7}, {-3, -1, 12, 0, 7}, {1, -1, 1, 3, 13}}; //det F = -710 printf("det A = %.16f\n", detTriang(3, A)); printf("det B = %.16f\n", detTriang(2, B)); printf("det C = %.16f\n", detTriang(4, C)); printf("det D = %.16f\n", detTriang(4, D)); printf("det E = %.16f\n", detTriang(4, E)); printf("det F = %.16f\n", detTriang(5, F)); return 0; }
Referências
- [1] CORMEN, T. H. et al. Algoritmos: teoria e prática. 3 ed. Rio de Janeiro: Elsevier, 2012.
- [2] DANTE, L. R. Matemática, Volume único. 1ª edição. São Paulo: Ática, 2005.
- [3] FRANCO, N. B. Cálculo Numérico. 1 ed. São Paulo: Pearson, 2006.
Nenhum comentário:
Postar um comentário