Algoritmo para decomposição LU

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\mathbb{R}^{n\times n}$, com todos os menores principais não singulares. Então $A$ pode ser decomposta de forma única em $$A=LU,$$ 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.

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=LU.$ 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} \displaystyle{\quad \sum_{j=1}^{n-1} (n-j)[1+2(n-j)] }\\ \displaystyle{\qquad \qquad=\sum_{j=1}^{n-1} n(1+2n) -j(1+4n) +2j^2}\\ \displaystyle{\qquad \qquad= n(1+2n)\sum_{j=1}^{n-1} 1 -(1+4n)\sum_{j=1}^{n-1} j + 2 \sum_{j=1}^{n-1}j^2}\\ \displaystyle{\qquad \qquad= n(1+2n)(n-1) - {1\over 2}n(n-1)(1+4n) + {2\over 3}(n-1)(n-1/2)n}\\ \displaystyle{\qquad \qquad= {2\over 3}n^3 + {\cal 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) é $${2\over 3}n^3 + {\cal O}(n^2) + {\cal O}(n^2) + {\cal O}(n^2) = {2\over 3}n^3 + {\cal 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. $$ \left[ \begin{array}{rrrr} 2 & 1 & 1 & 3 \\ -2 & -2 & 1 & 1 \\ 6 & 5 & 2 & 2 \\ -4 & -3 & 6 & 5 \end{array}\right] x = \left[ \begin{array}{r} -2 \\\ -3 \\\ 5 \\\ -9 \end{array} \right] \quad \mbox{e} \quad \left[ \begin{array}{rrrr} 2 & 1 & 1 & 3 \\ -2 & -2 & 1 & 1 \\ 6 & 5 & 2 & 2 \\ -4 & -3 & 6 & 5 \end{array}\right] x = \left[ \begin{array}{r} 0 \\\ -4 \\\ 0 \\\ -15 \end{array} \right]. $$

resposta

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

resposta

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?