Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como utilizar a decomposição LU para calcular o determinante de uma matriz?

+1 voto
471 visitas
perguntada Jun 13, 2017 em Matemática por Peng Yaohao (101 pontos)  
Compartilhe

1 Resposta

+1 voto
respondida Jun 13, 2017 por Peng Yaohao (101 pontos)  
editado Jun 30, 2017 por Peng Yaohao

A decomposição LU é uma fatoração que transforma uma matriz quadrada de ordem \(n\) no produto de duas outras matrizes, também \(nxn\), que são triangulares -- ou seja, matrizes que possuem todos os elementos acima (matrizes triangulares inferiores) ou abaixo da diagonal principal (matrizes triangulares superiores) iguais a zero. [Note que dizer que uma matriz é tanto triangular superior quanto inferior equivale a dizer que ela é uma matriz diagonal]. A decomposição LU tem esse nome justamente por isso: \(L\) vem de "lower triangular" e \(U\) vem de "upper triangular".

Basicamente, a decomposição LU reescreve uma matriz quadrada qualquer \(A_{nxn}\) na forma \(A_{nxn}=L_{nxn}U_{nxn}\), onde \(U\) é a forma escalonada de \(A\), obtida mediante aplicação sequencial de operações-linha, processo conhecido como Eliminação de Gauss, comumente usada para a resolução de sistemas lineares. São três as operações-linhas:

  1. \(L_i\leftrightarrow L_j; i,j\in\{1,2,...,n\}\) (permutação de linhas)
  2. \(L_i=k\cdot L_i;k\in\mathbb{R}; i,\in\{1,2,...,n\}\) (multiplicação por escalar)
  3. \(L_i=L_I-k\cdot L_j;k\in\mathbb{R};i,j\in\{1,2,...,n\}\) (soma com múltiplo de outra linha; \(a-b\) é simplesmente a soma \(a+(-b)\))

A decomposição LU pode ser utilizada como um método mais prático de se calcular o determinante de uma matriz quadrada \(A_{nxn}\) utilizando-se duas propriedades:

  • Para quaisquer A e B matrizes quadradas de mesma ordem, \(det(AB)=det(A)det(B)\);
  • Para qualquer matriz quadrada e triangular \(A_{nxn}, det(A)=\prod\limits_{i=1}^{n}{a_{ii}}\) (produtório dos elementos da diagonal principal)

Como por construção \(A=LU\), \(det(A)=det(L)det(U)\). Como \(L\) e \(U\) são matrizes triangulares, basta fazer o produto dos elementos de suas diagonais principais, um processo muito menos oneroso que a expansão por cofatores (método de Laplace). Como então encontrar \(L\) e \(U\)?

Sabe-se que a matriz \(U\) é a forma escalonada de \(A\), obtenível via operações-linha. Basicamente, busca-se zerar todos os elementos abaixo da diagonal principal de \(A\) utilizando uma das três operações-linha. Cada operação-linha está associada a uma matriz \(E_{nxn}\) chamada de matriz elementar, que é resultante de se aplicar a operação-linha na matriz identidade de ordem \(n\). É possível demonstrar que a eliminação de Gauss consiste num produto sequencial de \(m\) matrizes elementares a partir da matriz original \(A\), de modo que:

\((E_1E_2...E_m)\cdot A=U\).

É fácil observar que qualquer matriz elementar é inversível (basta aplicar a operação-linha inversa na matriz identidade). Assim, tem-se que:

\((E_m^{-1}...E_2^{-1}E_1^{-1})(E_1E_2...E_m)\cdot A=(E_m^{-1}...E_2^{-1}E_1^{-1})\cdot U\), donde resulta que \((E_m^{-1}...E_2^{-1}E_1^{-1})=L\)

A partir de uma matriz \(A\) qualquer, é fácil ver que a operação-linha 2 (multiplicação por escalar) não ajuda a obter os zeros abaixo da diagonal principal durante o escalonamento, pois o escalar em questão teria que ser \(0\), o que anularia a linha inteira. Assim, apenas as operações-linha 1 e 3 são usadas. A primeira é usada quando o elemento pivô já é zero, trocando a linha inteira por outra com valor não-nulo nessa coluna; e como ela não precisará ser zerada (pois já é zero), nenhuma operação adicional é exigida. Desse modo, a principal operação-linha é a terceira -- soma por múltiplo de outra linha: para se zerar o elemento \(a_{ij},i>j\), faz-se a operação-linha:

\(L_i=L_i-(\frac{a_{ij}}{a_{jj}})\cdot L_j\), onde \(a_{jj}\) é o elemento pivô.

Observa-se que a matriz elementar associada à operação-linha 3 é idêntica à matriz identidade, com exceção do elemento \(a_{ij}\), onde \(i\) e \(j\) são as linhas envolvidas na operação-linha, e que \(a_{ij}\) é justamente o fator \(\frac{a_{ij}}{a_{jj}}\) que aparece na operação-linha, com sinal trocado. Utiliza-se esse fato para construir a matriz triangular inferior \(L\) (o escalonamento definido aqui usa subtração, então o ``sinal trocado'' em relação à soma é exatamente \(\frac{a_{ij}}{a_{jj}}\)). Partindo-se da matriz identidade \(I\), a cada operação-linha número 3 executada em \(A\), atualiza-se os elementos \(i_{ij}\), substituindo-se os zeros pelos fatores \(\frac{a_{ij}}{a_{jj}}\) - como duas linhas \(i\) e \(j\) estão envolvidas, a diagonal principal de \(I\) permanece inalterada com valores todos iguais a \(1\). Após aplicar todas as operações-linha em \(A\) até chegar em \(U\), a matriz identidade atualizada será justamente a matriz \(L\)!

Para o cálculo do determinante, utilizam-se os fatos de que a operação-linha 3 não altera o valor do determinante, e a operação-linha 1 inverte seu sinal. Assim, a cada permutação de linhas, o sinal do determinante é trocado; ao fim do escalonamento, o determinante de \(A\) será simplesmente o produto dos elementos da diagonal principal de \(U\), já que a construção de \(L\) faz com que sua diagonal principal seja toda igual a \(1\).

Seguem implementações em Python dos algoritmos de expansão de Laplace e decomposição LU. A expansão de Laplace calcula o determinante de \(A\) a partir da expressão \(det(A)=\sum\limits_{k=1}^n{a_{ik}C_{ik}}\), onde \(C_{ij}=(-1)^{i+j}M_{ij}\) é o cofator, com \(M_{ij}\) sendo o menor - matriz resultante da exclusão da linha \(i\) e da coluna \(j\) de \(A\).

import numpy as np
import copy
import math

def laplace(A): ### Adaptado das notas de aula do Prof. Daniel Cajueiro (curso de métodos computacionais, Universidade de Brasília, 2017/1)
    n=np.shape(A)[0]
    if(n==1):
        det=A
    elif(n==2):
        det=A[0,0]*A[1,1]-A[0,1]*A[1,0]
    else:
        det=0
    for i in range(n):
        newMatrix=A[1:,:]
        newMatrix=np.delete(newMatrix,i,axis=1)
        det=det+math.pow(-1,1+i+1)*A[0,i]*(laplace(newMatrix))
    return det

def LU(x):
    dim=np.shape(x) # dimensões da matriz (linhas/colunas)
    det=1 # elemento neutro da multiplicação
    lower=np.identity(dim[0])  # começa-se pela matriz identidade para a matriz L
    for j in range(0,dim[1]):
        for i in range(j,dim[0]):  # eliminação de Gauss
            if(i == j and x[i,j]==0 and i!=dim[1]-1):  # permutação de linhas para deixari pivô diferente de zero
                p=-2  # contador provisório
                k=i
                while p!=k and k<dim[1]-1: # trocar linha pivô por outra que tenha elemento não nulo na coluna pivô
                    if x[k,j]!=0:
                        p=k # identifica primeira linha com coluna pivô não nula
                    else:
                        k+=1
                if p==-2: # caso não tenha (ou seja, coluna inteira nula), a matriz é singular
                    return("Determinante = 0 (Matriz singular)")
                else:
                    x[i,],x[p,]=x[p,].copy(),x[i,].copy() # permutação de linhas
                    det=(-1)*det # trocar sinal do determinante (pela propriedade)
            elif(i!=j): # com o pivô não-nulo:
                xt0=x[j,j] 
                xt1=x[i,j]
                ft=xt1/xt0
                x[i,:]=x[i,:]-x[j,:]*ft # operação-linha número 3
                lower[i,j]=ft # o fator que multiplica a outra linha é atribuída à matriz L - não se troca o sinal pois a operação já é de substração
        det=det*x[j,j] # após zerar todos os elementos abaixo do elemento pivô na respectiva coluna, agrupar o elemento da diagonal principal ao determinante - ao final de todas as colunas, o número resultante será o produtório de todos os elementos dessa diagonal, e portanto o determinante da matriz original (por construção, a matriz identidade de L será sempre 1)
    print "U=",x
    print "L=",lower
    print "LU = ",np.dot(lower,x), "= X"
    print "Determinante = ",det

A função "LU" irá retornar as matrizes \(L\) e \(U\), além do produto \(LU\) (que deve ser igual à matriz de entrada original) e do determinante. A complexidade desse algoritmo é \(\mathcal{O}(n^3)\), pois a obtenção de \(U\) (escalonamento) exige checar cada um dos \(n^2\) elementos da matriz, e para cada um deles, executar uma operação-linha. A atribuição dos elementos na matriz \(L\) e o produto agregado do termo de determinante têm custos fixos.

Em comparação com o método dos cofatores, esse método é bem menos dispensioso, pois a expansão de Laplace envolve calcular \(n\) determinantes de matrizes de ordem \(n-1\) até se atingir o caso base, de modo que para matrizes de ordem maior o crescimento do número de operações necessárias será exorbitante - por exemplo, uma matriz \(5x5\) exigiria calcular o determinante de 5 matrizes \(4x4\), cada uma das quais exigiria calcular 4 determinantes de matrizes \(3x3\), e assim por diante. Assim, o método dos cofatores possui complexidade \(\mathcal{O}(n!)\), significantemente maior que a decomposição LU.

Segue código para replicação:

import random

n=3 # n pode ser qualquer número inteiro positivo
x=np.empty([n,n],dtype=float ) # cria matriz nxn
for i in range(0,n):
    for j in range(0,n):
        x[i,j]=random.uniform(-10,10) # atribui números aleatórios de \(-10\) a \(10\) para cada elemento da matriz
LU(x) # aplica a função LU

Vamos agora comparar empiricamente a eficiência dos dois algoritmos, testando a velocidade de cada um deles para matrizes aleatórias de ordem \(nxn\), com \(n\) variando de \(1\) a \(11\):

import timeit

lim=11 # tamanho máximo da matriz
cont=np.empty([lim])
tempo=np.empty([lim,2])
for i in range(0,lim):
    cont[i]=i+1 # eixo x do gráfico
    test =  np.empty([i+1,i+1],dtype = float )
    for j in range(0,i+1):
        for k in range(0,i+1):
            test[j,k] = random.uniform(-100,100) # matrizes aleatórias de ordem 1 a 11
    start=timeit.default_timer()
    laplace(test)
    stop=timeit.default_timer()
    tempo[i,0]=stop-start # tempo de processamento força bruta (O(n!)), eixo y do gráfico azul
    start=timeit.default_timer()
    LU(test)
    stop=timeit.default_timer()
    tempo[i,1]=stop-start # tempo de processamento LU (O(n^3)), eixo y do gráfico vermelho

O algoritmo de força bruta demorou mais de \(4000\) segundos para calcular o determinante de uma matriz \(11x11\), enquanto que a mesma matriz demandou apenas \(0.0135\) segundos utilizando decomposição LU! Os resultados do teste para \(n\) variando de 1 a 11 estão sintetizados na figura abaixo:

### Adaptado das notas de aula do Prof. Daniel Cajueiro (curso de métodos computacionais, Universidade de Brasília, 2017/1)
import matplotlib.pyplot as plt
ax = plt.subplot(111)
ax.plot(cont,tempo[:,0],color='b')
ax.plot(cont,tempo[:,1],color='r')
ax.set_ylabel('Tempo de processamento (segundos)')
ax.set_xlabel('Ordem da matriz')
plt.show()

A imagem será apresentada aqui.

Valores maiores de \(n\) não foram testados por conveniência, mas é fácil observar que o tempo esperado para o cálculo do determinante de uma matriz \(12x12\) será \(12\) vezes maior que o tempo do cálculo de uma \(11x11\), algo em torno de \(48000\) segundos (mais de 13 horas!), uma \(13x13\) exigiria \(13\) vezes mais tempo, e assim por diante. Para uma matriz \(20x20\), demoraria aproximadamente 8 milhões de anos para calcular seu determinante usando a expansão por cofatores, enquanto que, usando a decomposição LU, o tempo de processamento foi de apenas \(0.3908\) segundo!

comentou Jun 30, 2017 por João Gabriel Souza (76 pontos)  
Ótimo exercício: apresenta de forma didática a aplicação do método LU, fazendo várias elucidações ao longo do exercício. A parte teórica do problema foi muito bem exposta e elucidativa. O exercício mostrou-se corretamente replicável, apresentando código comentado o que facilita o entendimento da replicação.

Um outro exemplo de decomposição L e U de matrizes pode ser encontrado no seguinte endereço: https://rosettacode.org/wiki/LU_decomposition#Python.

Apenas para fazer um adendo ao apresentado no exercício, é interessante observar que o método de eliminação de Gaus, utilizando no exemplo para a obtenção das matrizes L e U, é utilizado na resolução de problemas de Programação Linear. Um exemplo disso é o algorítimo SIMPLEX. O texto deixa claro a utilização da eliminação de Gaus para solução de sistemas lineares. Esse seria apenas um adendo ao exemplo exposto.
...