Fatorações LU e Cholesky

Transcript

1. Fatorações LU e Cholesky Prof. Paulo R. G. Bordoni UFRJ
2. Na scipy.linalg o método de eliminação de Gauss é o solve( ).
3. Um programa usando o solve( ):
4. E uma execução dele.
5. Nas próximas transparências mostraremos que o método de eliminação de Gauss pode ser descrito com operações matriciais. Teremos então a fatoração LU.
6. Ao final da sequência de multiplicações, como na eliminação, obteremos uma matriz triangular superior. As operações elementares usadas no algoritmo de eliminação de Gauss (com pivotamento) podem ser descritas por multiplicações matriciais !
7. | = 11 12 ⋯ 1 1 22 ⋯ 2 2 ⋱ ⋮ ⋮ | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋯ ⋮ Sim, da matriz aumentada ] obteremos a matriz ] abaixo. A letra de (superior em inglês), uma denominação tradicional.
8. Mostre como trocar duas linhas da | de posição, Loirinha! É mais fácil que trocar o pneu do meu carro !
9. É isso mesmo Loirinha. É só multiplicar pela matriz identidade com linhas trocadas. O produto da matriz identidade com as linhas r e s trocadas é a matriz com essas linhas trocadas.
10. Como sou objetivo, fiz um programa para conferir.
11. Uma execução do meu programa:
12. As operações elementares no processo de eliminação foram realizadas da seguinte forma: Linha 1 Linha 1 Linha 2 Linha 2 Linha 1 21 /11 × − Linha 3 Linha 3 Linha 1 31 /11 × − Linha n Linha n Linha 1 1 /11 × − Multiplicadores
13. 1 = − 1 = 1 0 ⋯ 0 0 1 ⋮ ⋮ 0 ⋯ ⋱ 0 0 1 − 0 0 ⋯ 0 21 0 ⋮ ⋮ 1 ⋮ 0 ⋱ 0 ⋯ 0 = 1 0 ⋯ 0 −21 1 ⋯ 0 ⋮ −1 0 ⋱ ⋮ ⋯ 1 , com 1 = Τ 1 11 , = 2, ⋯ , . Pois é, na 1ª etapa isto corresponde a efetuar o produto 1 ∙ = ( − 1 ) ∙ .
14. Este programa que mostra o produto da = − por A.
15. Seu programa é show, Mestra !
16. Exatamente, Cabelos de Fogo. Agora é só multiplicar essa matriz E0 ∗ Aum por E1 para obter o 2º passo do processo de eliminação. É o que farei !
17. Tomei o cuidado de mostrar passo-a-passo o processo de eliminação matricial.
18. Observem que (|)2 = (2 |)1 = (2 1 |)0 .
19. E, veja só, Loirinha: (|)4 = (|) onde é matriz triangular superior do final do processo de eliminação ! Portanto (|)4 = 4 (|)3 = = 4 3 (|)2 = 4 3 2 (|)1 = = 4 3 2 1 (|)0
20. Resumindo = 4 3 2 1 (|) . E portanto (4 3 2 1 )−1 = (|).
21. Mestre, todas as matrizes são triangular inferior, logo o produto 4 3 2 1 delas também será. E a inversa (4 3 2 1 )−1 também ! Lembrei-me que a inversa de um produto de matrizes é o produto das inversas na ordem trocada. Assim (4 3 2 1 )−1 = 1 −12 −13 −14 −1
22. E, algo surpreendente, a inversa de cada matriz de eliminação = − é a matriz + , obtida simplesmente trocando o sinal da : Essa é forte! Eu tenho que conferir!
23. Meu programa:
24. Confirmou, Surfista !
25. O que é verdadeiramente surpreendente Surfista, é que o produto = 1 −12 −13 −14 −1 é superposição das matrizes dos multiplicadores !!! Faça um programa para calcular a matriz Surfista. E, comprove que ela é uma matriz com 1′ na diagonal e que abaixo da diagonal estão os multiplicadores usados em cada etapa do processo de eliminação.
26. E que = , isto é, que = . Este resultado é conhecido como Teorema da fatoração LU. Resumindo tudo, acabamos de verificar que o processo de eliminação é realizado pela matriz triangular inferior = 1 −12 −13 −14 −1.
27. A paisagem da Scipy.linalg é magnífica; há muita coisa a ser conhecida!
28. Álgebra linear Apresentação Rotinas básicas Decomposições Funções matriciais Matrizes especiais Agora vamos explorar as fatorações (Decompositions) existentes na scipy.linalg.
29. A famosa fatoração LU.
30. Rotinas existentes na Referência da scipy.linalg para a fatoração LU.
31. Sim, de = , segue que det() = det det det(). Veremos que L só possui 1’s na diagonal e que P é ortogonal. Assim, det = ± det = ±11 22 ⋯ A fatoração = permite calcular det().
32. Mestre, que outras vantagens Maria Loirinha leva com essa tal de fatoração LU? Antes do Mestre te responder, jovem, proponho que discuta com seus colegas sobre essa cultura de levar vantagem em tudo. Grande Filósofo!
33. Respondendo tua pergunta, Loirinha: Tendo a fatoração = , resolver o sistema linear = é equivalente a resolver = . Então resolvemos primeiro o sistema = e depois o = . Note que esses dois são sistemas triangulares, facílimos de resolver! Conforme já vimos, o maior trabalho reside em fatorar (i.é, efetuar o processo de eliminação).
34. Os detalhes para usar a linalg.lu( ):
35. O que linalg.lu( ) retorna:
36. A rotina pode efetuar o pivotamento por linha (ou não). Ela retorna: • as três matrizes P, L e U • ou apenas duas, a PL que é a triangular inferior permutada e a U.
37. Conferindo a fatoração LU:
38. Um 1º exemplo em que não necessidade de trocas de linhas,
39. Já neste exemplo é evidente a necessidade de permutar a L0 com a L2 e depois a L2 com a L3. Atente para as matrizes A e P.
40. Uma outra forma de chamar a fatoração LU da SciPy. Como no 1º exemplo, não houve a necessidade de pivotamento.
41. Chamando a 2ª forma da fatoração LU da SciPy. Agora com pivotamento.
42. Esta é lu_factor( ) que efetua a fatoração . Ela recebe a matriz A e retorna uma matriz contendo a fatoração LU e um vetor P contendo os índices de pivoteamento.
43. Os parâmetros retornados. Esta rotina usa a LAPACK.
44. Esta rotina recebe a fatorada LU e o vetor piv fornecidos pela lu_factor( ). Recebe também o termo independente b e retorna a solução x.
45. Resolvendo um sistema linear Ax = b via fatoração LU com pivotamento.
46. Ela é boa para resolver sistemas Ax=b com muitos b’s, mas a mesma A.
47. Vejam, entrei com uma matriz 4x4 e 2 termos independentes. Poderiam ser muitos: 10, 50, etc.
48. Uma matriz = [ ] possui banda inferior k quando = 0 para − > . E possui banda superior k, quando = 0 para − > . = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × Eis uma matriz A com banda inferior 2 (em laranja) e banda superior 1 (em verde):
49. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ Se A é uma matriz 100x100 teríamos 10.000 elementos para armazenar. Nesse armazenamento de banda precisaremos apenas de 4x100-4 ∎ ∎ ∎ ∎ A ação para armazenar uma matriz de banda como essa é “esquecer os zeros e deixar cair”! Vejam só:
50. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ ∎ ∎ ∎ ∎ Mestre, prefiro pensar em “deitar” as diagonais.
51. Esta é rotina específica da scipy.linalg para resolver sistemas lineares cuja matriz é de banda.
52. Novamente, os detalhes sobre os parâmetros
53. Um programinha usando a solve_banded( )
54. Um outro programa usando a solve_banded( )
55. Uma execução do programa.
56. Um exemplo maior! Nesse caso A possui 72 = 100-28 zeros.
57. Quando é igual à sua transposta: = . Simétrica Quando satisfaz: ≠ 0 ⟹ = ෍ ,=1 > 0 Definida positiva Quando sua transposta é a inversa: = −1 ⟺ = = Ortogonal Continuando a visita ao zoológico matricial:
58. 1. + = + 2. ()= 3. = 4. ( )= 5. = Observações: • De 1 e 2 decorre que matrizes simétricas constituem um subespaço do espaço de todas as matrizes. • Atenção: a 5. não permite concluir que a matriz produto é uma matriz simétrica. Algumas propriedades da transposição:
59. De fato, o produto de matrizes simétricas não é uma matriz simétrica! E, além disso, o produto também não é comutativo. Vejam os contra-exemplos! 3 2 2 1 0 1 1 −2 = 2 −1 1 0 0 1 1 −2 3 2 2 1 = 2 1 −1 0
60. A matriz A = 1 0 0 0 −2 0 0 0 3 não é definida-positiva, pois se = temos: = 1 0 0 0 −2 0 0 0 3 = = −2 3 = 2 − 22 + 32 e escolhendo = 2 1 0 obtemos = 0.
61. A matriz = 5 −1 2 −1 2 2 2 2 3 é definida-positiva, pois para = temos: = − 2 + + 2 2 + 2 + 2, uma soma de quadrados, que só se anula para = = 0 0 0 . Surfista, confira essa afirmação da Mestra!
62. = sen cos 0 −cos sen 0 0 0 1 sen − cos 0 cos sen 0 0 0 1 = 1 0 0 0 1 0 0 0 1 Loirinha, agora é tua vez de conferir a afirmação da Mestra! A matriz = sen − cos 0 cos sen 0 0 0 1 é ortogonal pois :
63. Um zoológico bem complexo: A adjunta de = [ ], anotada AH, é definida por = [ത ]. É a conjugada da transposta. Adjunta Quando satisfaz: = . É a correspondente de uma matriz simétrica, com elementos complexos, Hermitiana, ou auto-adjunta Quando sua adjunta é a inversa: = −1 ⟺ = = Unitária
64. As matrizes ortogonais (e as unitárias, no caso complexo) desempenham um papel fundamental em Álgebra Linear Computacional. Aguardem!
65. Esta rotina é apropriada para matrizes hermitianas (simétricas no caso real) e positivo-definidas.
66. Detalhamento dos parâmetros.
67. O exemplo de armazenamento do Tutorial da scipy.linalg. Armazena-se apenas a diagonal e a banda superior (ou inferior). Em vez de 36 = 6 × 6 elementos, apenas 15 = 3 × 6 − 3.
68. Esta 1ª parte do programa recebe uma matriz A (formato de banda), um termo independente b e devolve a solução x usando a rotina para matrizes simétricas, positivo-definidas e de banda da scipy.linalg. O nome do programa é: sist_banda_sim_def_pos.py
69. Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
70. Uma execução do programa. Se a matriz A não for definida positiva, o programa acusa. Vejam a seguir.
71. O 4º elemento da diagonal principal é o responsável por tornar a matriz A não positiva definida. A linalg acusa o fato! ERRO
72. Uma forma de calcular a inversa de uma matriz A é resolvendo a equação matricial = usando solve( )
73. Uma brevíssima explicação do que é a fatoração de Cholesky.
74. Não deixe de ver no livro do Trefethen.
75. É essa a ideia Surfista. Além disso, vale a volta; veja na próxima transparência! Confirmando, se A é simétrica e positiva definida está garantida a existência de uma matriz triangular superior tal que = . É como se U fosse a raiz quadrada de A, Mestre, pois = .
76. Suponha que uma matriz não-singular A possa ser fatorada na forma = , com U triangular superior. Então A será simétrica e definida positiva. É difícil provar isto Mestra?
77. E também que A é definida positiva, isto é ≠ 0 ⟹ > 0. De fato assumindo ≠ 0 e = teremos ≠ 0, devido à não-singularidade de U. Então = () = = = 2 > 0. Bem, minha querida, de = segue que: • A é simétrica, pois = ()= ()= = • U é não-singular, pois det()2 = det() = det() ≠ 0.
78. As rotinas do SciPy para a fatoração de Cholesky.
79. A explicação para passagem de parâmetros.
80. Surfistas e Loirinhas, fiz este programa usando a fatoração de Cholesky.
81. Uma execução. Deixo a cho_fator() e a cho_solve() para vocês.
82. Esta rotina só efetua a fatoração. Deve ser usada com a cho_solve().
83. Fácil, fácil. Só troquei o nome da rotina e acrescentei o parâmetro de retorno Eh_inf.
84. Fácil nada, você errou. O apressado come cru! Repeti os dados da Mestra e não conferiu.
85. Errou porque não prestou atenção no aviso: Vou corrigir!
86. Foi só aplicar a rotina triu() na Cho. Vejam o resultado a seguir:
87. Como você disse Surfista, fácil, fácil!
88. Esta rotina resolve o sistema após fatorado pela rotina cho_factor().
89. Vou resolver o sistema linear abaixo. A matriz do sistema é a do exemplo anterior. 5 2 1 0 2 7 3 1 1 0 3 1 6 −1 −1 2 = 1 2 3 4
90. O Mestre chamou a Cho_fator() e depois a Cho_solve(). Percebam que ele não precisou zerar a parte inferior da matriz Cho. Além disso, b entrou como vetor-linha e foi transformado em vetor-coluna.
91. A execução, com conferência a resposta.
92. Esta é específica para matrizes de banda. Só fatora.
93. Esta resolve o sistema linear, após a fatoração pela cholesky_banded().
94. Façam um programa para resolver um sistema linear = no qual , além de ser simétrica e positiva-definida, também é de banda. Vai cair na prova Mestra?
95. Tchau, até a próxima aula!