Algoritmo para decomposição LU

Construção e análise do algoritmo para a decomposição LU

Vamos lá?

Na aula anterior vimos como o processo de escalonamento deu origem à decomposição LU, através de um exemplo. Nesta aula vamos escrever o algoritmo para a decomposição LU e analisá-lo.

Sempre é possível computar a decomposição LU de uma matriz quadrada?

A. Sim, o algoritmo apresentado faz exatamente isto.
B. Não, o algoritmo apresentado pode falhar em computar a decomposição LU.

Se a matriz não tiver decomposição LU, o que acontecerá ao aplicar o algoritmo?

A. O algoritmo identificará que a matriz não tem decomposição LU.
B. Uma decomposição errada será computada, dado que a matriz não cumpre as hipóteses do algoritmo.

Qual o custo da decomposição LU, em termos no número de operações de ponto flutuante?

A. É da ordem de $2 n^{3}$ operações.
B. É da ordem de $n^{2}$ operações.
C. É da ordem de $\frac{2}{3} n^{3}$ operações.
D. Tem custo total de $\frac{2}{3} n^{3} + 2 n^{2}$ operações.

O experimento numérico do final do vídeo...

A. refutou o cálculo do número de operações de ponto flutuante feito para o algoritmo.
B. validou o cálculo do número de operações de ponto flutuante feito para o algoritmo.
C. mostrou que o tempo de execução não é proporcional ao número de operações de ponto flutuante no algoritmo.
D. indicou que outros fatores também contribuem para o tempo de execução, além do número de operações de ponto flutuante do algoritmo.

Nem toda matriz tem decomposição LU. O teorema a seguir caracteriza exatamente quais matrizes a têm, em termo de seus menores principais. Um menor principal de ordem $j$ de $A$ é a submatriz de ordem $j$ , formada apenas pelas primeiras $j$ linhas e colunas, ou seja, a submatriz quadrada de ordem $j$ no canto superior esquerdo de $A$ .

Teorema: Seja $A \in R^{n \times n}$ , com todos os menores principais não singulares. Então $A$ pode ser decomposta de forma única em $A = L U,$ onde $L,$ de ordem $n,$ é triangular inferior com diagonal unitária e $U,$ de ordem $n,$ é triangular superior.

Talvez você fique com a impressão de que antes de tentar aplicar o algoritmo para computar a decomposição LU é necessário verificar se matriz atende às condições do teorema. Será? Observe o algoritmo para decomposição LU apresentado no vídeo.

$U \leftarrow A,$ $L \leftarrow I$
Para $j = 1, 2, \dots, n - 1$

Se $u_{j j} = 0,$ então $A$ não tem decomposição LU. FIM.
Para $i = j + 1, \dots, n$

$l_{i j} \leftarrow u_{i j} / u_{j j}$
$u_{i j} \leftarrow 0$
Para $k = j + 1, \dots, n$

$u_{i k} \leftarrow u_{i k} - l_{i j} u_{j k}$

Na prática, caso a matriz $A$ tenha decomposição LU, este algoritmo produzirá duas novas matrizes ao final de sua execução, as matrizes $L$ e $U,$ tais que $A = L U .$ Porém, caso a matriz $A$ não tenha decomposição LU, o algoritmo vai identificar isto por um pivô nulo (única possibilidade que inviabiliza a existência da decomposição) e encerrar a execução. Desta forma, muito mais fácil (e barato!) que tentar verificar se a matriz $A$ satisfaz ou não o teorema, é aplicar diretamente o algoritmo acima.

Na aula passada, vimos que um sistema linear qualquer pode ser resolvido em duas etapas: (1) computa-se a decomposição LU e (2) resolvem-se dois sistemas lineares triangulares. Já vimos em outra aula o algoritmo para a resolução de um sistema linear triangular superior e contamos o número de operações de ponto flutuante realizadas. Para saber o custo total da resolução de um sistema linear geral, resta computar quantas operações de ponto flutuante são consumidas no algoritmo acima.

No passo 2.2.3.1, duas operações são realizadas. Como esse passo é repetido para $k$ de $(j + 1)$ a $n,$ temos um total de $2 (n - j)$ operações ou flops. Mais uma operação é realizada no passo 2.2.1. Esse montante é repetido, por conta do loop do passo 2.2, de $(j + 1)$ a $n .$ Desta forma, o custo completo do passo 2.2 é $(n - j) [1 + 2 (n - j)] .$ Por fim, o loop do passo 2 itera de $j = 1$ até $j = (n - 1) .$ Aqui não podemos simplesmente multiplicar o número de repetições pelo custo interno do loop, uma vez que este custo interno varia com o índice $j .$ Sendo assim, temos que o custo completo é dado por $\begin{array}{l} \sum_{j = 1}^{n - 1} (n - j) [1 + 2 (n - j)] \\ = \sum_{j = 1}^{n - 1} n (1 + 2 n) - j (1 + 4 n) + 2 j^{2} \\ = n (1 + 2 n) \sum_{j = 1}^{n - 1} 1 - (1 + 4 n) \sum_{j = 1}^{n - 1} j + 2 \sum_{j = 1}^{n - 1} j^{2} \\ = n (1 + 2 n) (n - 1) - \frac{1}{2} n (n - 1) (1 + 4 n) + \frac{2}{3} (n - 1) (n - 1 / 2) n \\ = \frac{2}{3} n^{3} + O (n^{2}) . \end{array}$

Portanto, o custo para resolver o sistema linear original (soma do custo da decomposição LU com o custo dos dois sistemas triangulares) é $\frac{2}{3} n^{3} + O (n^{2}) + O (n^{2}) + O (n^{2}) = \frac{2}{3} n^{3} + O (n^{2}) .$ Ou seja, face ao custo da decomposição LU, o custo da resolução dos sistemas lineares triangulares é desprezível.

Para verificar se de fato o tempo computacional para o cálculo da decomposição é LU é diretamente proporcional ao número de operações de ponto flutuante, foi realizado um experimento numérico. Para matrizes aleatórias de diferentes dimensões, foi medido o tempo gasto na decomposição LU, computada pela rotina interna do Octave. Essa mediação de tempo foi repetida 100 vezes em cada experimento e o valor considerado foi a média de todas essas repetições. Os dados assim coletados estão no gráfico abaixo.

Usando o a técnica de ajuste de curvas por quadrados mínimos, descobriu-se que o tempo de cálculo da decomposição LU, em função da dimensão da matriz, é dado por $t (n) \approx (7.70 \cdot 10^{- 9}) n^{2.26} .$ Se o tempo fosse diretamente proporcional ao número de operações deveríamos ter encontrado o expoente $3$ ao invés de $2.26 .$ Isso ocorre porque há outros fatores envolvidos que impactam no tempo de execução, como a movimentação de dados entre a memória do computador e o processador e principalmente o emprego de algum grau de paralelismo (além, claro, da possibilidade do algoritmo implementado não ser exatamente o algoritmo exibido aqui). Mais a frente neste curso, há um projeto onde você terá a oportunidade de repetir este experimento e descobrir a curva que descreve o tempo de cálculo da decomposição LU em seu próprio ambiente computacional.

Para a próxima aula, veremos um exemplo computacional pequeno onde a resolução do sistema linear não sai como o esperado.

Referência

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

1. Compute a decomposição LU e utilize-a para resolver os sistemas lineares abaixo. $[\begin{array}{rrrr} 2 & 1 & 1 & 3 \\ - 2 & - 2 & 1 & 1 \\ 6 & 5 & 2 & 2 \\ - 4 & - 3 & 6 & 5 \end{array}] x = [\begin{array}{r} - 2 \\ - 3 \\ 5 \\ - 9 \end{array}] e [\begin{array}{rrrr} 2 & 1 & 1 & 3 \\ - 2 & - 2 & 1 & 1 \\ 6 & 5 & 2 & 2 \\ - 4 & - 3 & 6 & 5 \end{array}] x = [\begin{array}{r} 0 \\ - 4 \\ 0 \\ - 15 \end{array}] .$

A solução do primeiro sistema é $x = (2, - 1, 1, - 2)^{T}$ e a solução do segundo sistema é $x = (- 1, 2, - 3, 1)^{T} .$

2. Verifique que $A = [\begin{matrix} 3 & 1 & 0 \\ 6 & 2 & 1 \\ 1 & 1 & 1 \end{matrix}]$ é não singular. Mostre que $A$ não tem decomposição LU.

A matriz $A$ é não-singular uma vez que $det (A) = - 2,$ que pode ser computado facilmente apenas porque a matriz é de baixíssima ordem. Ao tentar computar a decomposição $L U,$ após a primeira iteração, o pivô para a segunda iteração é nulo e portanto a matriz não tem decomposição LU.

3. Considerando o seu computador, qual a ordem da maior matriz quadrada capaz de ser armazenada em memória (cada número em ponto flutuante ocupa 4 bytes)? Para simplificar, suponha que 100% da memória do seu computador está disponível para isso. Quantas operações de ponto flutuante seriam necessárias para computar a decomposição LU dessa matriz? Usando a fórmula para $t (n)$ encontrada no experimento numérico descrito nesta aula, quanto tempo seria necessário para computar a decomposição LU dessa matriz?

Confira Aprofunde Exercite Discuta