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

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

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

Por efetuar

Paulo Bordoni

November 07, 2014
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. Agora vamos usar os superpoderes das NumPy e SciPy, para

    dirigir o magnífico calhambeque da capa pelo caminho das pedras.
  2. 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!
  3. 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.
  4. 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]
  5. 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!
  6. Agora fatiaremos a coluna 2 de uma matriz 3x4: 0

    1 2 3 0 1 2 [ ∶, 2] E agora Surfista?
  7. 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]
  8. 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: ]
  9. 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 )
  10. 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.
  11. 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.
  12. 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.
  13. As honras da Alemanha à K. F. Gauss, por muitas

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

    descreveram o processo em linhas gerais! Só papo...
  21. 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:
  22. 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.
  23. | = 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!
  24. 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
  25. 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
  26. 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
  27. 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
  28. Fiz um programa para mostrar cada etapa do processo de

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

    Comparem com a anterior. Mais adiante, apresentaremos uma versão aperfeiçoada.
  32. 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.
  33. 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.
  34. 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.
  35. 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!
  36. 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.
  37. Os dados são os mesmos e o sistema é impossível.

    Entretanto o código em A oferece a solução abaixo!!! B A
  38. 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!
  39. 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!
  40. 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.
  41. 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!
  42. 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:
  43. 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!
  44. 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.
  45. 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!
  46. Agora posso fechar os parênteses. Vou retornar ao método da

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

    Isto significa que os multiplicadores perderão 10 dos seus 23 preciosos bits de precisão.
  49. 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.
  50. 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!
  51. 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.
  52. 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
  53. Este é um 1º exemplo de execução. Surfista, faça outros,

    semelhantes aos exemplos anteriores. 1 2 3
  54. 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.
  55. 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.
  56. 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:
  57. | = 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.
  58. 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.
  59. • 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!
  60. 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.
  61. 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 ⟷ .
  62. = 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.
  63. Ele permite conferir que o produto de permutações é comutativo.

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

    = . Observe, Surfista, que neste exemplo escolhi = .
  65. 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
  66. 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.
  67. 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
  68. 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:
  69. 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!
  70. = 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.
  71. −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
  72. 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!
  73. 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!
  74. 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.
  75. 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!
  76. 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?
  77. 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
  78. 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.