Decomposição LU
Como fatorar uma matriz em um produto de duas matrizes triangulares
Vamos lá?
Na aula anterior vimos como o processo de escalonamento pode converter um sistema linear e outro sistema equivalente, porém com a vantagem da matriz de coeficientes ser triangular superior. Nesta aula veremos que este processo pode ser feito em duas etapas, uma onde apenas a matriz é manipulada e outra, posterior, onde o vetor independente é utilizado.
- 1
- 2
- 3
- 4
No escalonamento da primeira coluna de $A,$ todas as operações realizadas foram da forma $\ell_i \leftarrow - ( m )\ell_1 + \ell_i,$ onde $\ell_1$ é a primeira linha de $A$ e $\ell_i$ é a linha a ser atualizada. Como o multiplicador $m$ foi computado?
O que acontece se em alguma etapa do processo o pivô for nulo?
A matriz $L$ obtida na decomposição LU...
Durante o processo de eliminação gaussiana para simplificar o sistema linear $Ax = b,$ operamos na matriz aumentada do sistema, $[A|b].$ Todas as operações realizadas alteram tanto as entradas da matriz $A$ quanto do vetor $b.$ Entretanto, se você olhar com atenção o exemplo da aula passada, vai reparar que apenas as entradas da matriz $A$ foram determinantes na definição de cada operação. O vetor $b,$ mesmo tendo sido alterado, em nenhum momento foi decisivo para a escolha da operação elementar a ser aplicada na matriz aumentada.
Esta constatação nos motiva a postergar a manipulação do vetor $b.$ Numa primeira fase, manipulamos apenas a matriz $A,$ triangularizando-a através de operações elementares, e numa fase seguinte, operamos no vetor $b,$ chegando assim a um sistema triangular superior, cuja resolução é simples. Uma vantagem imediata em separar o processo nessas duas etapas é que, caso queiramos resolver dois sistemas lineares com a mesma matriz $A$ mas vetores $b$ distintos, apenas a segunda etapa do processo precisa ser refeita, uma vez que a triangularização da matriz $A$ é comum aos dois sistemas.
Para tanto, as operações elementares realizadas em $A$ precisam ser registradas para posteriormente também serem aplicadas ao vetor $b.$
Cada operação realizada no processo de eliminação é da forma $$ \ell_i \leftarrow - ( m ) \ell_j + \ell_i, \quad i = j+1,\ldots,n, $$ onde $\ell_j$ é a linha do pivô, $\ell_i$ é a linha onde se quer eliminar o elemento da coluna do pivô e $n$ é a dimensão da matriz. O multiplicador $m$ é razão entre o elemento a ser eliminado e o pivô. Essa operação, pode ser representada pelo produto à esquerda por uma matriz especialmente montada como $$ \begin{array}{cc} \begin{array}{ccccccc} \hphantom{1} & \hphantom{\ddots} & j & \hphantom{\ddots} & i & \hphantom{\ddots} & \hphantom{1}\\ \hphantom{1} & \hphantom{\ddots} & \downarrow & \hphantom{\ddots} & \downarrow & \hphantom{\ddots} & \hphantom{1}\\ \end{array} & \\ \left[ \begin{array}{ccccccc} 1 \\ & \ddots \\ & & 1 \\ & & & \ddots \\ & & -m& & 1 \\ & & & & & \ddots \\ & & & & & & 1 \end{array} \right] & \begin{array}{l} \vphantom{1}\\ \vphantom{\ddots} \\ \leftarrow j \\ \vphantom{1} \\ \leftarrow i \\ \vphantom{\ddots} \\ \vphantom{1} \end{array} \end{array} $$
Essa matriz é tão especial que sua inversa é simples. Se o produto por esta matriz tem como efeito somar um múltiplo da linha pivô à linha $i,$ a matriz inversa deve desfazer essa operação. Ou seja, a inversa deve subtrair o mesmo múltiplo. Desta forma, para obter a inversa basta trocar o sinal do elemento na posição $(i,j),$ de $-m$ para $m.$ Além disso, o produto acumulado por matrizes como essa produz uma matriz que "sobrepõe" os todos os multiplicadores a uma matriz identidade. Desta forma, podemos decompor ou fatorar a matriz $A$ como o produto de duas matrizes $L$ e $U,$ onde $L$ é uma matriz triangular inferior com diagonal unitária, composta pelos multiplicadores computados durante o processo (apenas a divisão entre os elemento a ser eliminado e o pivô), e $U$ é uma matriz triangular superior. Essa decomposição (ou fatoração) é conhecida como decomposição LU.
Mas por que escrever $A=LU$ facilita a resolução do sistema linear original? Observe que se $$b = Ax = (LU)x = L (Ux).$$ Se denotarmos $Ux$ por $y,$ teremos que $Ax=b$ se torna $Ly = b,$ onde $Ux=y$, ou seja, temos dois sistemas triangulares encadeados para resolver. A solução do primeiro sistema, o vetor $y,$ é o termo independente para o segundo sistema. De fato, o vetor $y$ é o resultado de aplicar ao vetor $b$ todas as operações elementares computadas previamente no escalonamento de $A.$
Um ponto que não foi discutido no exemplo do vídeo é o que acontece se o pivô for, em algum momento do processo de escalonamento, nulo. Se isto acontecer, não há como usar a linha do pivô para eliminar as entradas abaixo da diagonal e portanto a matriz $A$ não tem decomposição $LU.$
Se do ponto de vista teórico, apenas o caso do pivô nulo é um problema, do ponto de vista numérico devemos esperar complicações quando o pivô for muito pequeno também. Entretanto, antes de entrar nesse assunto, na próxima aula vamos formalizar o algoritmo para o cálculo da decomposição LU e descobrir qual o custo computacional envolvido.
Essa aula foi baseada no excelente livro Trefethen e David Bau (1997).
Referência
Lloyd N. Trefethen e David Bau, III. Numerical Linear Algebra. SIAM, 1997.
1. Se $ L = \begin{bmatrix} 1 & 0 & 0\\ 1/2 & 1 & 0 \\ 1/3 & 1/2 & 1 \end{bmatrix}$ e $U= \begin{bmatrix} 2 & 4 & -3 \\ 0 & 3 & 1 \\ 0 & 0 & 3 \end{bmatrix} $ são os fatores da decomposição LU de uma matriz $A,$ resolva o sistema linear $Ax=b,$ para $b=(18, 20, 17.5)^T.$
Para resolver $Ax=b,$ usando os fatores $L$ e $U,$ resolve-se primeiro $Ly=b,$ cuja solução é $y=(18, 11, 6)^T,$ depois resolve-se $Ux=y,$ cuja solução é $x=(6, 3, 2)^T.$ De forma alguma deve-se formar a matriz $A$ pelo produto de $L$ por $U,$ para então resolver $Ax=b.$
2. Suponha que $A$ tem decomposição LU. Discuta como utilizar os fatores $L$ e $U$ para calcular (a) o determinante de $A,$ e (b) a matriz inversa, $A^{-1}.$