Upgrade to Pro — share decks privately, control downloads, hide ads and more …

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

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

Por efetuar.

Avatar for Paulo Bordoni

Paulo Bordoni

June 02, 2015
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. O título deste conjunto de transparências poderia ser “Pra não

    dizer que não falei de flores” Canção escrita e interpretada por Geraldo Vandré. Ficou em segundo lugar no Festival Internacional da Canção (1968 – os Anos de Chumbo) . Teve sua execução proibida durante anos, após tornar-se um hino de resistência do movimento civil e estudantil que fazia oposição à ditadura militar brasileira.
  2. Agora vamos usar os superpoderes das NumPy e SciPy, para

    dirigir o magnífico calhambeque da capa pelo caminho das pedras. Falar de flores... Nenhum aluno nosso sairá deste curso de Cálculo Numérico sem saber “Eliminação de Gauss” e “fatoração LU”.
  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. Esta notação é idêntica à do Matlab. [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] E agora Surfista?
  8. 0 1 2 3 0 1 2 Se A é

    uma matriz m x n, então: • [ , : ] fornece a linha i de A, a partir de j; • [ : , ] fornece a coluna j de A, de i para baixo, inclusive i. [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] : k `` = ´´ [ j, k )
  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. Confira pelas listas dos índices das linhas e das colunas.

    As cores mostram a ordem. 0 1 2 3 0 1 2 4
  13. Confira pelas listas dos índices das linhas e das colunas.

    Atente para a ordem das colunas. 0 1 2 3 0 1 2 4 3 4 5
  14. 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.
  15. 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.
  16. As honras da Alemanha à K. F. Gauss, por muitas

    outras contribuições à Matemática. E as nossas também, Mestre!
  17. 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.
  18. 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.
  19. 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 ~ .
  20. 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.
  21. 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.
  22. 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.
  23. Sim! Como é feita, de fato, a eliminação? Mestres, vocês

    descreveram o processo em linhas gerais! Só papo...
  24. 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:
  25. 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.
  26. | = 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!
  27. 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
  28. 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
  29. 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
  30. 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
  31. Linha 1 Linha 1 Linha 2 Linha 2 Linha 1

    21 /11 × − Linha n Linha n Linha 1 1 /11 × − Sim Sherlock, na 1ª etapa: • Efetuamos − 1 divisões, para calcular os multiplicadores, • Em cada linha − 1 multiplicações e subtrações • Como são − 1 linhas: • ( − 1)2 multiplicações • ( − 1)2 subtrações
  32. 1 1 + 2 2 … = {1 1 +

    2 2 … } − (1 /11 ){ 11 1 + 21 2 + ⋯ } 1 = 0 Multiplicações Atenção Surfista, como já sabemos que 1 = 0 não precisamos efetuar essa conta! Discordo, são n multiplicações e n subtrações! Subtrações
  33. 11 1 + 12 2 + ⋯ + 1 =

    1 22 2 + ⋯ + 2 = 2 ⋯ 2 2 + ⋯ + = Depois repetimos o processo para o subsistema sem a 1ª linha e assim por diante! Assim é só somar as operações em todas as n-1 etapas. Veja a conta a seguir, Surfista!
  34. No processo todo, para uma matriz = [ ] de

    ordem teremos: • ( − 1) + ( − 2) + ⋯ + 1 divisões, • ( − 1)2+( − 2)2+ ⋯ + 12 multiplicações, • ( − 1)2+( − 2)2+ ⋯ + 12 subtrações. Em resumo, o número de operações envolvidos na redução da matriz = [ ] de ordem n a uma matriz triangular superior é: • (23 − 32 + ) 6 multiplicações e subtrações, • (2 − 3 + 2) 2 divisões.
  35. Para essa conclusão usamos que: • 1 + 2 +

    ⋯ + = ( + 1)/2 • 12 + 22 + ⋯ + 2 = + 1 (2 + 1)/6 Vou tentar provar essas afirmações, Mestra!
  36. Sim! Como é feito, de fato, o processo de eliminação?

    Torno a repetir, vocês são bons de papo... E no computador?
  37. Fiz um programa para mostrar cada etapa do processo de

    eliminação, seguindo a ideia do Surfista. Vejam o resultado.
  38. 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
  39. 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
  40. Aos meus pés uma versão sem o loop em i.

    Comparem com a anterior. Mais adiante, apresentaremos uma versão aperfeiçoada.
  41. 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.
  42. 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.
  43. 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.
  44. 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!
  45. 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.
  46. Os dados são os mesmos e o sistema é impossível.

    Entretanto o código em A oferece a solução abaixo!!! B A
  47. 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!
  48. 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!
  49. 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.
  50. 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!
  51. 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:
  52. 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!
  53. 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.
  54. 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!
  55. Agora posso fechar os parênteses. Vou retornar ao método da

    eliminação de Gauss e focar todos os holofotes no pivot.
  56. 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 .
  57. Multiplicadores 210 × 21 210 × 31 210 × 1

    Isto significa que os multiplicadores perderão 10 dos seus 23 preciosos bits de precisão.
  58. 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.
  59. 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!
  60. 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.
  61. 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
  62. Este é um 1º exemplo de execução. Surfista, faça outros,

    semelhantes aos exemplos anteriores. 1 2 3
  63. 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.
  64. 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.
  65. 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:
  66. | = 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.
  67. 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.
  68. • 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!
  69. 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.
  70. 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 ⟷ .
  71. = 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.
  72. Ele permite conferir que o produto de permutações é comutativo.

    Quero lembrar que, em geral, o produto matricial não é comutativo.
  73. Ele permite, ainda, mostrar que as permutações são idempotentes: 2

    = . Observe, Surfista, que neste exemplo escolhi = .
  74. 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
  75. 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.
  76. 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
  77. 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:
  78. 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!
  79. = 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.
  80. −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
  81. 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!
  82. 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!
  83. 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.
  84. 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!
  85. 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?
  86. 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
  87. 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.