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:
- \(L_i\leftrightarrow L_j; i,j\in\{1,2,...,n\}\) (permutação de linhas)
- \(L_i=k\cdot L_i;k\in\mathbb{R}; i,\in\{1,2,...,n\}\) (multiplicação por escalar)
- \(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()
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!