LU com pivoteamento

Para reduzir erros numéricos quando resolvemos um sistema linear em aritmética de precisa finita, devemos evitar trabalhar com pivôs pequenos durante o processo de escalonamento. Isso vale tanto para a eliminação de Gauss, quanto para a decomposição LU. Um exemplo de como a magnitude relativa do pivô tem influência na solução numérica do sistema linear foi visto na aula passada.

A forma tradicional de se fazer isso é através da estratégia de pivoteamento parcial. Relembre que o algoritmo para decomposição LU sem pivoteamento é:

Neste algoritmo, quando iniciamos o passo $j,$ ou seja, quando formos eliminar os elementos abaixo da diagonal da coluna $j,$ o pivô é $u_{jj}$ (nesse ponto, $U$ é a matriz temporária, obtida no processo de escalonamento das colunas anteriores). A estratégia de pivoteamento parcial consiste em trocar a linha do pivô $j,$ por alguma das linhas abaixo dela, digamos a linha $k,$ determinada de forma que \begin{equation}\label{pivocond} |u_{kj}| =\max_{j\le i \le n} |u_{ij}|. \end{equation} Feita a troca das linhas (caso tenha sido preciso), o algoritmo segue como antes. Para registrar essa troca de linhas, precisamos de uma nova matriz, denominada matriz de permutação $P.$ Essa matriz nada mais é que a matriz identidade, sobre a qual se realizam as mesmas operações de troca de linhas que tenham sido realizadas em $U.$

O algoritmo para a decomposição LU com pivoteamento parcial fica assim

Vamos detalhar alguns passos do algoritmo acima. No passo 2.1 determinamos a linha $k,$ abaixo da linha $j$, que contem o maior elemento em módulo abaixo da diagonal. Se $k\neq j,$ então será necessário fazer o pivoteamento. Isto é feito nos passos de 2.2.1 a 2.2.3. Por simplicidade, estamos usando aqui a notação própria do Octave, a saber $A(j,:)$ representa a linha $j$ da matriz $A$ e $A(j,p:q)$ representa a linha $j$ de $A,$ apenas entre as colunas $p$ e $q.$ Os passos 2.2.1 e 2.2.2 são apenas as trocas de linhas de $U$ e $P$ como já havíamos comentado. Mas o passo 2.2.3 merece atenção. Nesse passo, estamos trocando as linhas de $j$ e $k$ da matriz $L,$ mas apenas da primeira coluna a até a coluna $j-1,$ pois apenas essas foram computadas nesse momento da execução do algoritmo. Essa troca é necessária pois cada elemento da matriz $L$ armazena um multiplicador que estava relacionado a duas linhas no processo de escalonamento (a linha pivô e a linha onde a operação de eliminação foi aplicada). Caso as linhas troquem de ordem por conta do pivoteamento, os multiplicadores também precisam ser trocados, pois do contrário essas associação entre multiplicadores e linhas ficariam incorretas.

No passo 2.3, após o pivoteamento, verificamos se o pivô é nulo. Se este for o caso, isto significa que todos os elementos abaixo da diagonal também são nulos (por quê?). Logo, não há nada mais a ser eliminado e podemos passar para o processo de eliminação dos elementos na próxima coluna.

Perceba que agora não há mais a possibilidade do algoritmo interromper sua execução antes do fim e, como consequência, toda matriz tem decomposição LU com pivoteamento parcial. Isso é diferente do que acontecia com decomposição LU sem pivoteamento, que pode existir ou não.

Ao final da execução do algoritmo acima, temos três matrizes $P,$ $L$ e $U$, que satisfazem $$ PA=LU. $$ Para resolver um sistema linear $Ax=b,$ premultiplicamos ambos os lados dessa equação por $P,$ assim o produto $PA$ pode ser substituído por $LU.$ Chamando $Ux=y$ temos $$ \left\{ \begin{array}{l} L y = Pb, \\ Ux = y. \end{array} \right. $$ Ou seja, o sistema original é reescrito como dois sistemas lineares, um triangular inferior e outro triangular superior.

Antes de concluir, um último comentário sobre a implementação computacional deste algoritmo. Caso a matriz original $A$ não seja mais necessária após sua decomposição nos fatores $L,$ $U$ e $P,$ podemos utilizar o espaço de memória que abrigava $A$ para armazenar os fatores $L$ e $U.$ Como essas duas matrizes são triangulares, uma inferior e outra superior, a porção triangular inferior de $A$ abrigará $L,$ enquanto que a porção triangular superior abrigará $U.$ Claro que falta espaço para armazenar os elementos da diagonal de $L$ e $U.$ Porém, a diagonal de $L$ é constante e igual a $1,$ não precisando ser armazenada. Com essa estratégia, não é necessário espaço adicional para os fatores $L$ e $U$ e o algoritmo acima fica até um pouco mais compacto (ainda mais usando operações vetoriais do Octave!).

Na próxima aula, veremos um exemplo prático, computado passo a passo, onde a decomposição LU com pivoteamento parcial é computada e e empregada na resolução do sistema linear.

Referência

Lloyd N. Trefethen e David Bau, III. Numerical Linear Algebra. SIAM, 1997.

1. Se $ L = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1/5 & 1 & 0 & 0 \\ -2/5\hphantom{-} & 3/4 & 1 & 0 \\ 4/5 & 2/4 & -1/2\hphantom{-} & 1 \end{array}\right],$ $U= \left[ \begin{array}{rrrr} 5 & -2 & 2 & 1 \\ 0 & 2 & 1 & -5 \\ 0 & 0 & 1 & 6 \\ 0 & 0 & 0 & -3 \end{array}\right]$ e $P =\begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} $ são os fatores da decomposição LU com pivoteamento de uma matriz $A,$ resolva o sistema linear $Ax=b$, para $b=(-5, 5, -7, -4.5)^T.$

2. Aplique o algoritmo da decomposição LU com pivoteamento parcial à matriz $$ A = \left[ \begin{array}{rrrr} -3.0 & 3.0 & 1.50 & -9.0 \\ 6.0 & 2.0 & -1.00 & 0.0 \\ 1.5 & -1.5 & 0.25 & 4.5 \\ 1.0 & -1.0 & 0.50 & 6.0 \end{array} \right]. $$ Usando os fatores da decomposição, resolva o sistema linear $Ax=b$, para $b = (9, 8, -6.5, -8)^T.$ Quando encontrar a solução, verifique se ela está correta, substituindo-a no sistema linear.

3. Reescreva o passo 2.4.3 do segundo algoritmo apresentado usando a notação compacta do Octave. Você deve escrever esse passo sem usar um loop.