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

Análise detalhada do processo de eliminação de Gauss e das razões IEEE 754 para pivotamento. Depois de explicar fatiamento, construímos programas vetorizados para resolver sistemas lineares, sem e com pivotamento. Em seguida mostramos como obter a fatoração LU a partir do processo de eliminação.

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, : ]
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 3 por 2:
16. Multipliquei uma sub-matriz por 0.5:
17. Fiz o produto direto de uma subcoluna por uma submatriz:
18. O produto direto de entre duas matrizes compatíveis é realizado em NumPy/SciPy através do produto delas como ndarrays. Portanto elas precisam ser convertidas antes, usando o modificador . A (esse A é de array): nome_da_matriz . A
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 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 k Linha 1 Antiga linha k 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 k Linha k 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 k Nova linha k 11 ≠ 0
36. Brilhante esta sua observação Surfista! Aliás adorei este filme.
37. Linha 1 Linha 1 A NumPy é mais poderosa ainda! Vejam: [1: , : ] = 1, : − ([1,1: ]/11 ) × [1, : ] Linha 2 Linha 3 Linha n − × Linha 2 Linha 3 Linha n 21 /11 31 /11 1 /11 Linha 1 Linha 1 Linha 1 11 ≠ 0
38. O 11 é intitulado “pivot”. Veremos adiante que, como os pivots de basquete, ele será o maior da coluna. Essas operações exigem que o pivot seja não-nulo: 11 ≠ 0. Caso 11 = 0 interrompemos o processo. 1 /11 31 /11 Multiplicadores 21 /11 11 ≠ 0
39. O final do processo é o mesmo: uma matriz triangular superior. Repetimos apenas para evidenciar que não precisamos envolver as incógnitas no processo – só os vetores linha da matriz aumentada. | = 11 12 ⋯ 1 1 22 ⋯ 2 2 ⋱ ⋮ ⋮ Observem que todos os elementos da matriz triangular aumentada são, eventualmente, diferentes dos iniciais.
40. Esta é a 1ª versão do código para o método de Gauss:
41. Dois exemplos onde o processo é interrompido.
42. Uma interrupção forçada. Abro parênteses para notícias cinzentas!
43. Um outro motivo é oriundo da maneira como números reais são tratados no computador. A representação IEEE 754 de ponto flutuante, no padrão single, será ilustrada na próxima transparência. Evitar divisão por zero não é a única razão para efetuar trocas de linhas!
44. 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:
45. 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!
46. 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.
47. 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!
48. Agora posso fechar os parênteses. Vou retornar ao método da eliminação de Gauss e focar os holofotes no pivot.
49. 1 /11 31 /11 Multiplicadores 21 /11 Assumam que o valor do pivot 11 é 1/1000. 11 ≅ 2−10 Então ao dividirmos por 11 estaremos multiplicando por 1.000 (≅ 1024 = 210) os elementos 21 , 31 , … , 1 .
50. Multiplicadores 210 × 21 210 × 31 210 × 1 Isto significa que os multiplicadores perderão 10 dos seus 23 preciosos bits de precisão.
51. 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.
52. Uma forma de minimizar a propagação de erro consiste em localizar o 1 de maior valor absoluto e trocar a linha k 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!
53. Usem esta função para buscar a posição do pivot
54. Os detalhes da argmax( ) para matrizes.
55. Os detalhes da argmax( ) para ndarrays.
56. Vejam elas funcionando:
57. O código para pivotamento por coluna só muda onde marquei.
58. 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 + ⋯ + =
59. | = 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.
60. 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.
61. • 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!
62. Começaremos pelas trocas, Loirinha!
63. 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.
64. 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 ⟷ .
65. = 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.
66. Este programa faz o que está descrito no seu início. Ele também mostra que o produto de permutações também é uma permutação.
67. Mestre, corrija o erro que assinalei, não é QPA, mas sim PQA ! He, he, he! Eu fiz como pegadinha, só para lembrar que = , uma consequência de serem matrizes de permutação.
68. Surfista, você já reparou que as matrizes de permutação são matrizes ortogonais? É verdade, Mestre! No seu exemplo, as colunas de PQ são a base canônica de ℝ5 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
69. 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.
70. Meu programa também permite conferir que −1 = .
71. 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.
72. O programa a seguir mostra isto:
73. As operações elementares 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
74. 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:
75. Um programa que mostra o produto da 1 = − 1 por A.
76. 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!
77. O programa que fiz para conferir a afirmação da Mestra!
78. = 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.
79. −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
80. Eu fiz este outro programa para conferir essa 2ª afirmação Mestra!
81. Vejam o resultado!
82. Tchau, até a próxima aula!