Da eleiminação de Gauss à fatoração PLU

Iniciamos recordando o operador de fatiamento, depois descrevemos o processo de eliminação de Gauss e como programá-lo com linalg. Mostramos como decidir se um sistema apresenta uma, nenhuma ou infinitas soluções. Apontamos a necessidade da realização do pivotamento face à aritmética finita da IEEE 754. Em seguida mostramos como o processo pode ser realizado através de multiplicações matriciais. Fechamos com a fatoração LU e a PLU.

Transcript

1. LNCC Da eliminação de Gauss à fatoração LU Prof. Paulo R. G. Bordoni UFRJ
2. Agora vamos usar os superpoderes das NumPy e SciPy, para dirigir o magnífico calhambeque da capa pelo caminho das pedras.
3. 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!
4. 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.
5. 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. [1, : ] 0 1 2 3 0 1 2 3 [ ∶, 2]
6. 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!
7. Agora fatiaremos a coluna 2 de uma matriz 3x4: 0 1 2 3 0 1 2 [ ∶, 2]
8. 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]
9. 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: ]
10. 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]
11. 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.
12. Agora vejamos como pegar submatrizes, como a abaixo: 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
13. Passei uma matriz M e peguei uma submatriz. 0 1 2 3 0 1 2 3 4 5 [1: 3,2: 5]
14. Agora vou passar manteiga nas fatias de pão, Mestra.
15. Multipliquei um pedaço da coluna 1 por 1.5:
16. Multipliquei uma sub-matriz por 0.5:
17. Fiz o produto direto de uma subcoluna por uma submatriz:
18. 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.
19. 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.
20. Fui buscar estas informações na Wikipedia!
21. As honras da Alemanha à K. F. Gauss, por muitas outras contribuições à Matemática. E as nossas também, Mestre!
22. 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.
23. 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.
24. 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 ~ .
25. 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.
26. 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.
27. 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. Alguns textos falam em escalonamento.
28. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
29. 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:
30. 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.
31. | = 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
32. 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
33. 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
34. 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
35. Mestre, com 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
36. Brilhante esta sua observação Surfista! Aliás adorei este filme.
37. Fiz um programa para mostrar cada etapa do processo de eliminação, seguindo a ideia do Surfista. Vejam o resultado.
38. Este é o código do programa. Marquei a observação do Surfista.
39. 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
40. 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
41. Aos meus pés uma versão sem o loop em i. Comparem com a anterior. À frente, apresentaremos uma versão aperfeiçoada.
42. Neste 1º exemplo o sistema possui uma única solução.
43. 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.
44. 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.
45. 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.
46. Uma interrupção forçada. Abro parênteses para notícias cinzentas!
47. 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”
48. A B Veja a seguir um excelente motivo para usar código confiável, como scipy.linalg, BLAS e LAPACK.
49. A B Os dados são os mesmos e o sistema é impossível. Entretanto o código em A oferece a solução abaixo!!!
50. 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!
51. 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!
52. 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.
53. O Mestre apresentará, a seguir, um “flash” da representação IEEE 754 de ponto flutuante. Mais adiante, no curso, faremos um mergulho em profundidade. Evitar divisão por zero não é a única razão para efetuar trocas de linhas!
54. 1 bit para o sinal 4 bytes = 32 bits 8 bits para o expoente 23 bits para a fração O padrão single na representação IEEE 754 possui 32 bits, distribuídos da seguinte forma:
55. 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!
56. Acreditem: 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 afirmar que a precisão de um single é de 23 bits.
57. Multiplicar por 2 desloca esse erro p casas binárias para a esquerda. Perdem-se p bits na precisão! Impressionante Sherlock, nunca tinha ouvido falar disto!
58. Agora posso fechar os parênteses. Vou retornar ao método da eliminação de Gauss e focar todos os holofotes no pivot.
59. 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 .
60. Multiplicadores 210 × 21 210 × 31 210 × 1 Isto significa que os multiplicadores perderão 10 dos seus 23 preciosos bits de precisão.
61. 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.
62. 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 pivoteamento parcial por coluna. Na realidade esta é, apenas, “a ponta do iceberg”. Analisaremos tudo isto em detalhes, mais adiante!
63. Usem esta função para buscar a posição do pivot
64. Os detalhes da argmax( ) para matrizes.
65. Os detalhes da argmax( ) para ndarrays.
66. Vejam elas funcionando:
67. Mestre, eis como trocar duas linhas de uma matriz:
68. 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.
69. 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
70. Este é um 1º exemplo de execução. Surfista, faça outros, semelhantes aos exemplos anteriores.
71. Mestre, 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.
72. Valeu, Mestre!
73. Acabamos de ver a programação do método da eliminação, de Gauss, para obtenção da solução = 1 2 ⋯ de um sistema de equações lineares: 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + =
74. | = 11 12 ⋯ 1 ⋮ 1 22 ⋯ 2 ⋮ 2 ⋱ ⋮ ⋮ ⋮ ⋮ | = 11 12 ⋯ 1 ⋮ 1 21 22 ⋯ 2 ⋮ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ ⋮ ⋮ ⋮ Efetuamos operações elementares sobre as linhas da matriz aumentada | até transformar A numa matriz triangular superior U.
75. 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.
76. • P é uma matriz de permutações (p/. trocas de linhas) • L é uma matriz triangular inferior (L de lower) • U é uma matriz triangular superior (U de upper) A fatoração PLU é o resultado final dessas multiplicações. Todas essas três operações elementares podem ser descritas através de multiplicações matriciais!
77. Começaremos pelas trocas, Loirinha!
78. 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.
79. 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 ⟷ .
80. = 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.
81. Este programa faz o que marquei aqui no início.
82. Ele permite conferir que o produto de permutações é comutativo. Quero lembrar que, em geral, o produto matricial não é comutativo.
83. Ele permite mostrar, também, que a inversa da permutação ↔ , é a permutação ↔ .
84. Ele permite, ainda, mostrar que (↔ )2 = , i. é, as permutações são idempotentes: neste exemplo = .
85. 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
86. 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.
87. A multiplicação da linha k de uma matriz A por um número α também é descrita por um produto matricial. A matriz M obtida a partir da identidade I substituindo o 1 da posição (i, i) por α faz isto.
88. O programa a seguir mostra isto:
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 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 1 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!
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 mostraremos a fatoração LU sem pivotamento.
100. Mostramos passo a passo as matrizes Mult, M e b. Ao final a solução x e o produto LU.
101. As mudanças no código da eliminação de Gauss com pivotamento, para mostrar a fatoração PLU são pequenas. Comparem com o código anterior.
102. Observe que o produto PLU, que já traz as trocas de linhas incluídas.
103. Tchau, até a próxima aula!