Da eliminação de Gauss à fatoração LU

Por efetuar

Transcript

1. LNCC Da eliminação de Gauss à fatoração LU Prof. Paulo R. G. Bordoni UFRJ
2. O título deste conjunto de transparências poderia ser “Para não dizer que não falei de flores”
3. Agora vamos usar os superpoderes das NumPy e SciPy, para dirigir o magnífico calhambeque da capa pelo caminho das pedras.
4. Antes, precisaremos aprender a obter fatias de ndarrays e a operar com elas. Na NumPy, a faca para fatiar arrays é o operador : . É a ferramenta que ela oferece para cortar o pão!
5. Representarei uma matriz por uma grade, como abaixo. Por que só estou interessada nos índices. Isto basta para explicar o conceito de fatiamento. 0 1 2 3 0 1 2 3 5 Lembrem-se, em NumPy, a indexação de arrays começa por zero e não por um.
6. 0 1 2 3 0 1 2 3 5 Se A é uma matriz m x n, então: • [ , : ] fornece a linha i de A; • [ ∶ , ] fornece a coluna j de A. Esta notação é idêntica à do Matlab. [1, : ] 0 1 2 3 0 1 2 3 [ ∶, 2]
7. Nossa 1ª fatiada será a linha 1 de uma matriz 4x5: 0 1 2 3 0 1 2 3 5 [1, : ] Esta é outra opção, Mestra!
8. Agora fatiaremos a coluna 2 de uma matriz 3x4: 0 1 2 3 0 1 2 [ ∶, 2] E agora Surfista?
9. 0 1 2 3 0 1 2 Se A é uma matriz m x n, então: • [ , : ] fornece a linha i de A, de j em diante; • [ : , ] fornece a coluna j de A, de i para baixo. [1,2: ] 0 1 2 3 0 1 2 3 [ 1: , 0]
10. Agora a fatia é uma parte da linha 1 de A: da posição 2 até o final. Faça o código para a outra figura Surfista. 0 1 2 3 0 1 2 [1,2: ]
11. No fatiamento devemos pensar em intervalos fechados à esquerda, mas abertos à direita: • [ , : ] fornece a linha i de A, iniciando na coluna j e acabando na coluna k-1. • Pense em j:k como o intervalo [j, k). 0 1 2 3 0 1 2 3 4 5 [1,1: 4] : k `` = ´´ [ j, k )
12. 0 1 2 3 0 1 2 3 4 5 [1,1: 4] Uma parte da linha 1, aquela correspondente aos índices [1,4) das colunas.
13. Agora vejamos como pegar submatrizes, como a abaixo: 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
14. Passei uma matriz M e peguei uma submatriz. 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
15. Agora vou passar manteiga nas fatias de pão, Mestra.
16. Multipliquei um pedaço da coluna 1 por 1.5:
17. Multipliquei uma sub-matriz por 0.5:
18. Fiz o produto direto de uma subcoluna por uma submatriz:
19. Portanto, antes, elas precisam ser convertidas, usando o modificador . A (esse A é de array): nome_da_matriz . A O produto direto de entre duas matrizes compatíveis é realizado em NumPy/SciPy através do produto delas como ndarrays.
20. Carl Friedrich Gauss 30/04/1777 – 23/02/1855 Passaremos a analisar o famoso Método de Eliminação, de Gauss. Na verdade ele não foi criado por Gauss. Ele já aparece no cap. 8 do livro chinês “Os nove capítulos sobre a arte matemática”, composto por gerações de mestres ao longo dos séculos X a II A.C.
21. Fui buscar estas informações na Wikipedia!
22. As honras da Alemanha à K. F. Gauss, por muitas outras contribuições à Matemática. E as nossas também, Mestre!
23. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Vamos agora descrever como resolver um sistema linear como o abaixo, pelo método de eliminação. Descreveremos o processo, primeiramente, em linhas gerais, usando as equações.
24. Em cada etapa utilizaremos operações elementares sobre as linhas. Elas são três: Somar ou subtrair linhas Multiplicar uma linha por um número real Trocar duas linhas de posição Nosso objetivo será atingido após várias etapas consecutivas.
25. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Na 1ª etapa, eliminaremos os termos envolvendo 1 nas equações após a 1ª, efetuando operações elementares nas linhas. 11 1 + 12 2 + ⋯ + 1 = 1 22 2 + ⋯ + 2 = 2 ⋯ 2 2 + ⋯ + = Os coeficientes e os termos independentes mudam após essas operações – por isto o ~ .
26. 11 1 + 12 2 + ⋯ + 1 = 1 22 2 + ⋯ + 2 = 2 ⋯ 2 2 + ⋯ + = Repetimos o processo para o subsistema sem a 1ª linha: Um sistema com n-1 equações a n-1 incógnitas.
27. Repetindo n-1 vezes o processo chegamos a um sistema triangular superior. 11 1 + 12 2 + ⋯ + 1 = 1 = 22 2 + ⋯ + 2 = 2 ⋯ Todos os coeficientes e termos independentes são diferentes dos originais, exceto, talvez, na 1ª linha.
28. O último passo consiste em resolver o sistema triangular. Algo extremamente simples. A parte trabalhosa do processo consiste, portanto, em eliminarmos incógnitas nas equações. Daí o nome: processo de eliminação. Os textos de Álgebra linear falam em escalonamento.
29. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
30. 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋯ 1 2 ⋮ = 1 2 ⋮ = Pois é... Agora vamos refazer tudo usando a forma matricial do sistema linear e a NumPy:
31. Acrescentando o vetor b à última coluna da matriz A, obtemos uma matriz referida na literatura como “matriz aumentada”. | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ A matriz aumentada contém todos os dados do problema.
32. | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ Como veremos, o processo de eliminação consiste em efetuar operações elementares sobre as linhas da matriz aumentada. Linha 1 Linha 1 Linha n Ahá! Por isso o fatiamento!
33. Se 11 ≠ 0, as operações elementares são descritas por: 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
34. 1 1 + 2 2 … = {1 1 + 2 2 … } − (1 /11 ){ 11 1 + 21 2 + ⋯ } 1 = 0 Nova linha i Linha 1 Antiga linha i Multiplic. Observem que os multiplicadores foram escolhidos de forma a anular todos os novos 1 abaixo do 11 . 1 /11 31 /11 Multiplicadores 21 /11 11 ≠ 0
35. As operações elementares são operações vetoriais: a multiplicação por fator de escala e a adição de vetores! Linha i Linha i Linha 1 −1 /11 × + Número real Vetor Vetor Vetor 11 ≠ 0
36. Mestre, com o fatiamento da a NumPy essa operação com linhas é feita num golpe só, como um Samurai! [, : ] = , : − (1 /11 ) [1, : ] Multiplic. Linha 1 Antiga linha i Nova linha i 11 ≠ 0
37. Brilhante esta sua observação Surfista! Aliás adorei este filme.
38. Fiz um programa para mostrar cada etapa do processo de eliminação, seguindo a ideia do Surfista. Vejam o resultado.
39. Este é o código do programa. Marquei a observação do Surfista.
40. Linha 1 Linha 1 A NumPy é mais poderosa ainda! Não precisamos do loop em i. Vejam todo o 1º passo de uma só vez! [1: , 1: ] = 1: , 1: − ([1: , 1]/11 ) × [1,1: ] Linha 2 Linha 3 Linha n − × Linha 2 Linha 3 Linha n 21 /11 31 /11 ,1 /11 Linha 1 11 ≠ 0
41. O 11 é intitulado “pivot”. Conforme veremos adiante, ele será o maior da coluna. Como os pivots de basquete. Essas operações exigem que o 11 seja não-nulo: 11 ≠ 0. Caso 11 = 0 interrompemos o processo. 1 /11 31 /11 Multiplicadores 21 /11 11 ≠ 0
42. Aos meus pés uma versão sem o loop em i. Comparem com a anterior. Mais adiante, apresentaremos uma versão aperfeiçoada.
43. Neste 1º exemplo o sistema possui uma única solução.
44. Neste 2º exemplo o sistema não possui uma única solução. Na etapa 2, a 4ª equação é 0 = −10. . Portanto o sistema é impossível.
45. Neste 3º exemplo o sistema não possui uma única solução. Na etapa 1, a 3ª equação é 0. = 0. . Portanto o sistema possui uma infinidade de soluções.
46. Estabelecendo uma primeira conclusão: quando possível, o final do processo é o mesmo, um sistema triangular superior para resolver. | = 11 12 ⋯ 1 1 22 ⋯ 2 2 ⋱ ⋮ ⋮ Senão, temos que examinar, também, como termo independente evolui no processo de eliminação.
47. Uma interrupção forçada. Abro parênteses para notícias cinzentas!
48. Apontarei claramente, agora, a diferença entre um curso de Cálculo Numérico e algo tipo: “Me engana que eu gosto, você finge que ensina e eu que aprendo”. Um detalhe aparentemente sem importância!
49. A B Veja a seguir um excelente motivo para usar código confiável, como scipy.linalg, BLAS e LAPACK. Assinalei a única diferença.
50. Os dados são os mesmos e o sistema é impossível. Entretanto o código em A oferece a solução abaixo!!! B A
51. Mestre, eu nunca teria posto esse detalhe em meu programa! Fiquei assustada. Como eu faço? Shht! Acabei de perder totalmente a confiança no meu taco!
52. Há mais de 40 anos o Rio entrava em choque com o desabamento do Elevado Paulo de Frontin, no Rio Comprido, Zona Norte do Rio. A tragédia causou a morte de 28 pessoas, deixou 30 feridas, e destruiu 17 carros. Foto: CPDoc JB Pois é jovens. O momento de aprender é este. Não permitam que os enganem! Imaginem se a consequência do erro tivesse sido esta!
53. O que causou o desastre no Paulo de Frontin eu não sei, mas o erro apresentado no código A deve-se ao desconhecimento da representação de números no computador, estabelecida pelo padrão IEEE 754 de ponto flutuante, internacionalmente adotado para computação numérica.
54. O Mestre recordará, a seguir, como um “flash” a representação IEEE 754 de ponto flutuante. Evitar divisão por zero não é a única razão para efetuar trocas de linhas!
55. 1 bit para o sinal 4 bytes = 32 bits 8 bits para o expoente 23 bits para a fração Lembrem-se, o padrão IEEE 754 single tem seus 32 bits, distribuídos na forma:
56. No single, a fração é constituída por 23 bits: 1 2 ⋯ 23 com = 0 1 . 23 bits para a fração E o valor da fração é dado por: = 1 2−1 + 2 2−2 + ⋯ + 23 2−23, sendo o restante dos bits desprezado!
57. Lembrem-se: A maioria esmagadora dos números reais na representação IEEE 754 de ponto flutuante, no padrão single, possui um erro na 24ª casa binária: 2−24. O Mestre acabou de reafirmar que a precisão de um single é de 23 bits.
58. Multiplicar por 2 desloca esse erro p casas binárias para a esquerda. Perdem-se p bits na precisão! Impressionante Sherlock, nunca tinha pensado nisto!
59. Agora posso fechar os parênteses. Vou retornar ao método da eliminação de Gauss e focar todos os holofotes no pivot.
60. 1 /11 31 /11 Multiplicadores 21 /11 Assumam que o valor do pivot 11 é 1/1000 = 0,001. 11 ≅ 2−10 Então ao dividirmos por 11 estaremos multiplicando por 1.000 (≅ 1024 = 210) os elementos 21 , 31 , … , 1 .
61. Multiplicadores 210 × 21 210 × 31 210 × 1 Isto significa que os multiplicadores perderão 10 dos seus 23 preciosos bits de precisão.
62. Que horror Mestra. É por isso que os pivots tem que ser (jogadores) grandes! Sim minha Loirinha querida! Faça como eu, despreze pivots pequenos, eles estragam o jogo, no final. Sempre vá em busca dos maiores possíveis.
63. Uma forma de minimizar a propagação de erro consiste em localizar o 1 de maior valor absoluto e trocar a linha i onde ele está com a linha 1. Esta é a estratégia conhecida na literatura como pivotamento parcial por coluna. Na realidade esta é, apenas, “a ponta do iceberg”. Analisaremos tudo isto em detalhes, mais adiante!
64. Usem esta função para buscar a posição do pivot
65. Os detalhes da argmax( ) para matrizes.
66. Os detalhes da argmax( ) para ndarrays.
67. Vejam elas funcionando:
68. Mestre, eis como trocar duas linhas de uma matriz:
69. Agora apresentaremos o código do método de eliminação de Gauss com pivotamento parcial por coluna. Precisamos localizar a posição do pivot na coluna, em seguida realizar a troca das linhas, para depois efetuar a eliminação.
70. O código só muda onde marquei. Em: 1. -> Busco o elemento de módulo máximo da coluna. 2. -> Faço a troca das linhas. 1 2
71. Este é um 1º exemplo de execução. Surfista, faça outros, semelhantes aos exemplos anteriores. 1 2 3
72. Mestre, como posso saber quais foram as trocas de linhas efetuadas? Muito simples, Loirinha. Basta criar um vetor de inteiros, trocas = [0,1, … , − 1] e anotar as trocas nele. Veja o código. O exemplo está na próxima transparência.
73. Valeu, Mestre! 1 3 2
74. Nas próximas 31 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 PLU.
75. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Lembrem-se, escrevemos todos os dados do sistema linear na matriz aumentada:
76. | = 11 12 ⋯ 1 ⋮ 1 22 ⋯ 2 ⋮ 2 ⋱ ⋮ ⋮ ⋮ ⋮ | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Depois efetuamos operações elementares sobre as linhas da matriz aumentada | até transformar A numa matriz triangular superior U.
77. As operações elementares usadas foram: Somar ou subtrair linhas Multiplicar uma linha por um número real Trocar duas linhas de posição E elas permitiram transformar o sistema original A x = b num sistema triangular superior U x = β, com a mesma solução x.
78. • P é uma matriz de permutações (p/. trocas de linhas) • L é uma matriz triangular inferior (L de lower) • U é a matriz triangular superior ao final da eliminação (U de upper) A fatoração PLU é o resultado final dessas multiplicações. Como o Mestre anunciou todas essas três operações elementares podem ser descritas através de multiplicações matriciais!
79. Começaremos pelas trocas, Loirinha!
80. Indicaremos por ↔ a matriz obtida permutando as linhas i e j da identidade I Observem que ↔ é uma matriz simétrica, isto é ↔ = ⟷ . É fácil observar que o produto ↔ é a matriz A com suas linhas i e j trocadas. O programa a seguir mostrar isto.
81. Neste programa você entra com uma matriz A pelo teclado e com os índices i e j das linhas de A a serem trocadas. Ele retorna a matriz ⟷ e o produto ⟷ .
82. = 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 Neste exemplo, mostro a identidade I com as linhas 1 e 2 e as 3 e 5 trocadas. Aliás, se P é uma matriz obtida da identidade I através de permutações de linhas, então o produto P A é a matriz A com as mesmas linhas trocadas.
83. Este programa faz o que marquei aqui no início.
84. Ele permite conferir que o produto de permutações é comutativo. Quero lembrar que, em geral, o produto matricial não é comutativo.
85. Ele permite mostrar, também, que a inversa da permutação ↔ , é a permutação ↔ .
86. Ele permite, ainda, mostrar que as permutações são idempotentes: 2 = . Observe, Surfista, que neste exemplo escolhi = .
87. Surfista, você já reparou que as matrizes de permutação são matrizes ortogonais? É verdade, Mestre! As colunas de uma permutação P são a base canônica de ℝ fora de ordem. = 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
88. Sim coleguinha. E da ortogonalidade segue que, se P é uma matriz de permutações então = −1. E como elas são simétricas temos −1 = = . Aliás isto é totalmente intuitivo: a inversa −1 destroca as linhas da identidade I que foram trocadas por P.
89. 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
90. 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 Pois é, isto corresponde a efetuar o produto 1 ∙ , Onde 1 = − 1 é a matriz de eliminação:
91. Um programa que mostra o produto da 1 = − 1 por A.
92. E, o mais surpreendente, é que a inversa dessa matriz de eliminação 1 = − 1 é a matriz + 1 , obtida simplesmente trocando o sinal da 1 : 1 −1 = ( − 1 )−1 = 1 0 ⋯ 0 21 1 ⋯ 0 ⋮ 1 ⋮ 0 ⋱ ⋮ ⋯ 1 = + 1 Essa é forte! Eu tenho que conferir!
93. O programa que fiz para conferir a afirmação da Mestra!
94. = 3 2 1 = 1 0 0 0 0 1 0 0 0 0 0 0 1 0 9 1 1 0 0 0 0 1 0 ⋮ 0 0 7 4 1 0 0 1 1 0 0 0 3 1 0 ⋮ −2 5 0 0 1 0 0 1 = = 1 0 0 0 3 1 0 0 19 188 7 67 1 0 9 1 No método de eliminação usamos um produto, de matrizes eliminadoras = 3 2 1 Como no exemplo abaixo.
95. −1 = 1 0 0 0 −3 1 0 ⋮ 2 −5 0 0 1 0 0 1 1 0 0 0 0 1 0 ⋮ 0 0 −7 −4 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 −9 1 = = 1 0 0 0 −3 1 1 0 2 −5 −7 −4 1 0 −9 1 O mais surpreendente é a inversa −1 = (3 2 1 )−1= 1 −12 −1 3 −1 desse produto E. Simplesmente é a superposição das com sinal trocado
96. Eu fiz este programa para efetuar a eliminação, passo a passo!
97. Vejam o resultado! A matriz A é atualizada em cada etapa.
98. Na realidade, ao programar, NÃO efetuamos multiplicações matriciais. É o que mostraremos nos dois programas a seguir. Então por que motivo mostraram tudo isto? Medite um pouco sobre tua pergunta, jovem Surfista imediatista!
99. Primeiro um programa para a fatoração LU sem pivotamento.
100. Ele mostra, passo a passo, as matrizes Mult, M e b. Ao final a solução x e o produto LU. Em cada etapa as matrizes Mult e M e o vetor b são atualizados!
101. As mudanças no código da eliminação de Gauss com pivotamento, para mostrar a fatoração PLU. Comparem com o código anterior.
102. Observe, Surfista, que o produto PLU já traz as trocas de linhas incluídas. Em cada etapa as matrizes M e L e o vetor b são atualizados!
103. Quando temos a fatoração = , resolver um sistema linear = é equivalente a resolver o sistema = , usando apenas a fatoração LU. Então Mestra, basta efetuar as trocas das linhas no termo independente b e usar a fatoração LU sem pivotamento na matriz A já com as trocas?
104. Sim, ao invés de resolver 1 1 −2 0 1 −5 2 −1 3 1 3 1 1 3 −2 2 0 1 2 3 = 3 −2 0 0 Resolvemos 3 3 1 −2 1 −5 2 −1 1 1 1 1 3 −2 2 0 0 1 2 3 = 0 −2 0 3
105. Foi o que eu fiz aqui. Usei a matriz M já com as linhas trocadas e troquei, também, as linhas do termo independente b. Depois resolvi o sistema sem pivotamento.
106. Tchau, até a próxima aula!