Cálculo de Determinantes via Triangularização

Por
|

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!).

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

CORMEN, T. H. et al. Algoritmos: teoria e prática. 3 ed. Rio de Janeiro: Elsevier, 2012.

DANTE, L. R. Matemática, Volume único. 1ª edição. São Paulo: Ática, 2005.

FRANCO, N. B. Cálculo Numérico. 1 ed. São Paulo: Pearson, 2006.

Siga o blog por e-mail

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