O determinante de uma matriz é uma função que, dada a matriz, retorna um número real. Na estatística, a própria fórmula da distribução normal multivariada apresenta o determinante da matriz de covariâncias.
O computador realiza a operação para nós, sem nos preocuparmos como está sendo calculado. Mas você já parou para pensar em, além de como é calculado, como é implementado o algoritmo? Vamos fazer isso nesse post, vamos implementar o cálculo de determinante pelo Teorema de Laplace que utiliza o método do cofator. Para isso, apresentarei duas formas de implementá-lo: o método recursivo e o método iterativo.
Se você já viu algo sobre recursividade em programação, sabe que ela pode ser um vilão para a memória do computador. Dessa forma, transformar o algoritmo recursivo em um algoritmo iterativo pode ser uma boa ideia. Nesse post faremos o algoritmo recursivo e, no próximo, o algoritmo iterativo.
A primeira coisa que vamos fazer é implementar a função que calcula o determinante de matrizes de ordem 2 e 3, que pode ser feito de forma simples e direta. É o que fazemos no seguinte código:
# Determinante de matrizes de ordem 2 e 3 ---------------
calc_det2 = function(A){
res = A[1,1]*A[2,2] - A[1,2]*A[2,1]
return(res)
}
calc_det3 = function(A){
res = A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2]
res = res - (A[1,3]*A[2,2]*A[3,1] + A[1,1]*A[2,3]*A[3,2] + A[1,2]*A[2,1]*A[3,3])
return(res)
}
Para matrizes de ordem maior que 3 o cálculo não é direto. Utilizamos, por exemplo, o teorema de Laplace, que diz que o determinante de uma matriz será dado pela soma do produto dos elementos de uma linha ou coluna (a sua escolha) pelos seus respectivos cofatores. Obviamente, esse método serve para as matrizes de ordem 3 e 2, porém o cálculo é mais complicado e o método direto é preferível.
Para calcular o cofator de um elemento, precisamos do determinante da matriz reduzida. De fato, $C(a_{ij}) = (-1)^{(i+j)} \times |D_{ij}|$, em que $a_{ij}$ é o elemento da linha i e coluna j; e $D_{ij}$ é a matriz obtida retirando-se a linha i e a coluna j (que estou chamando de matriz reduzida). Logo, necessitamos do determinante de uma outra matriz para calcular o determinante da matriz original. Por exemplo, se queremos calcular o determinante de uma matriz 4x4 e escolhemos a primeira linha para aplicar o método, precisamos calcular o determinante de quatro matrizes reduzidas 3x3 primeiro (uma para o cofator de cada elemento da linha escolhida).
O que isso está nos dizendo? Que se temos uma matriz n x n e queremos calcular seu determinante pelo método do cofator, precisamos calcular o determinante de todas as matrizes reduzidas de ordem 3 primeiro para depois calcular as de ordem maior, até chegar na matriz original.
Na verdade, poderíamos reduzir as matrizes até a ordem 1 no método, mas já sabemos calcular o determinante de matrizes de ordem 3 diretamente, então poupamos cálculo e também memória do computador.
Vamos começar nosso algoritmo. Declaramos nossa função:
calc_det = function(A){
...
}
Agora basta implementar. Pense comigo, se a matriz passada pelo parâmetro é de ordem 2 ou 3, já sabemos calcular, correto? Então, basta verificar a ordem da matriz e chamar uma das duas funções.
if(ncol(A)==3){
return(calc_det3(A));
}else if(ncol(A)==2){
return(calc_det2(A));
}
Agora, se a ordem da matriz é maior que 3, aplicamos o método do cofator. Podemos escolher qualquer linha ou coluna para aplicar o método e normalmente utilizamos a que contiver o maior número de zeros para facilitar os cálculos. Não faremos essa verificação aqui. Para simplificar o problema, vamos escolher sempre a primeira linha da matriz (inclusive das matrizes reduzidas, caso necessário) para realizar o método.
Além disso, precisamos de uma variável para armazenar o resultado do determinante e da matriz atual (matriz reduzida). Vejamos:
L1 = A[1,]
det_ = 0
auxM = A[-1,]
Até aqui estamos apenas preparando o terreno. Note que auxM é a matriz passada na função sem a primeira linha, como queríamos. Agora, temos que retirar as colunas, o que será feito de forma iterativa, isto é, para cada elemento da linha, retiramos sua respectiva coluna. Vejamos:
for(i in 1:length(L1)){
auxM2 = auxM[,-i]
auxdet = calc_det(auxM2)
cofat = (-1)^(1+i) * auxdet
det_ = det_ + L1[i]*cofat
}
Vamos por partes, das linhas internas ao for:
- Linha 1: calcula a matriz reduzida retirando a i-ésima coluna da matriz auxM.
- Linha 2: calcula o determinante da matriz reduzida.
- Linha 3: calcula o cofator associado ao elemento atual da linha escolhida da matriz.
- Linha 4: soma o cálculo do cofator ao resultado final.
Agora é a parte interessante. Quando o código vai da linha 2 para a linha 3? Apenas quando calculamos o determinante da matriz reduzida. Se você reparou, estamos chamando a função calc_det dentro dela. Chamamos isso de recursão em programação. Antes de avançar para a linha 3, a linha 2 deve ser finalizada, mas isso acontecerá apenas quando a função retornar o determinante da matriz reduzida.
Por exemplo, se estamos trabalhando com uma matriz 4x4, as matrizes reduzidas serão de ordem 3 e a função retornará o determinante dessas matrizes em uma etapa (pois calculamos diretamente). Mas se estamos trabalhando com uma matriz 5x5, as matrizes reduzidas serão de ordem 4 e, para calcular seu determinante, necessitamos das matrizes reduzidas de ordem 3. Logo, para avançar da linha 2 para a linha 3, vários códigos são executados recursivamente. Quando obtemos os determinantes das matrizes reduzidas de ordem 3, conseguimos calcular os determinantes das matrizes reduzidas de ordem 4 para, então, conseguir calcular o determinante da matriz original de ordem 5. Parece confuso, mas não é. Caso necessite, uma boa estudada sobre recursão ajuda.
Ao final da execução, devemos retornar o resultado:
return(aux)
O código completo da função está abaixo:
calc_det = function(A){
if(ncol(A)==3){
return(calc_det3(A));
}else if(ncol(A)==2){
return(calc_det2(A));
}
L1 = A[1,]
det_ = 0
auxM = A[-1,]
for(i in 1:length(L1)){
auxM2 = auxM[,-i]
auxdet = calc_det(auxM2)
cofat = (-1)^(1+i) * auxdet
det_ = det_ + L1[i]*cofat
}
return(det_)
}
Vamos testar nosso código. Primeiro com uma matriz de ordem 4:
D = matrix(c(1,2,3,4,2,1,3,4,2,1,9,7,4,5,1,4),ncol=4,byrow=T)
print(D)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 2 1 3 4
[3,] 2 1 9 7
[4,] 4 5 1 4
det(D)
[1] 72
calc_det(D)
[1] 72
A função det é a função implementada no R e estamos utilizando para verificar se nossa função calc_det está retornando o resultado correto. Podemos ver que nossa função retornou o resultado esperado.
Vamos testar com uma matriz de ordem 5:
E = matrix(c(3,9,1,8,1,7,4,9,2,1,2,3,
4,2,1,3,4,2,1,9,7,4,5,1,4),ncol=5,byrow=T)
print(E)
[,1] [,2] [,3] [,4] [,5]
[1,] 3 9 1 8 1
[2,] 7 4 9 2 1
[3,] 2 3 4 2 1
[4,] 3 4 2 1 9
[5,] 7 4 5 1 4
det(E)
[1] -629
calc_det(E)
[1] -629
Agora, temos nossa própria função para calcular o determinante de uma matriz via método do cofator (ou teorema de Laplace)!
No próximo post, faremos a versão iterativa do nosso algoritmo e vamos comparar a eficiência das função.
Espero que tenha gostado da aula.
Até a próxima aula!
