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

Álgebras linear computacional II

Paulo Bordoni
September 19, 2013

Álgebras linear computacional II

Recordamos os tipos de matrizes particulares usuais (simétricas, defin. posit., e de banda), apresentamos a utilização das rotinas básicas da sipy.linalg e suas fatorações LU e de Cholesky. Depois detalhamos a eliminação de Gauss e fatoração LU sem pivotamento e a necessidade do pivotamento para evitar a propagação do erro presente na aritmética finita da IEEE 754.

Paulo Bordoni

September 19, 2013
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. Pois é Surfista, a história mostra que, desde 1947, a

    comunidade de Computação Científica, tem trabalhado arduamente no desenvolvimento da Álgebra Linear Computacional. Em particular na resolução de sistemas lineares. Mestre, achei o caminho das pedras. Muito mais que isto, percebi claramente a importância das bibliotecas BLAS e LAPACK !
  2. Mestre, no tocante a resolução de sistemas lineares, proponho adotarmos

    a tática de operacionalizar primeiro. Depois vamos analisar os aspectos teóricos. Não gosto muito dessa ideia, mas como acredito na sua experiência didática, concordo!
  3. Voltando então às rotinas básicas do Tutorial, foquemos nossa atenção

    no tópico “Resolvendo sistemas lineares”. Rotinas básicas Achando a inversa Resolvendo sistemas lineares Achando o determinante Calculando normas Resolvendo problemas lineares de mínimos quadrados e pseudo-inversas Inversa generalizada
  4. Surfista, a orientação do Tutorial da scipy.linalg é usar a

    solve( ) [e não a inversa, inv( )] para obter a solução de um sistema linear.
  5. O que é estabilidade numérica, Mestre? Se é mais rápido

    é porquê envolve menos operações, mas como contá-las?
  6. O Surfista está certo, Mestra. Eis aí uma razão para

    não operacionalizarmos às cegas. Resmungão... Loirinha querida, mais adiante estudaremos erros, propagação de erros e estabilidade numérica.
  7. Recortei as rotinas básicas na Referência do scipy.linalg para resolução

    de sistemas lineares. Além da solve( ), existem rotinas apropriadas para matrizes com característica especiais.
  8. Aliás existe toda uma fauna de matrizes = [ ]

    cujas definições devemos recordar: Quando todos os seus elementos fora da diagonal são nulos: ≠ ⟹ = 0 Diagonal Quando todos os elementos acima da diagonal são nulos: < ⟹ = 0 Triangular inferior Quando todos os elementos abaixo da diagonal são nulos: > ⟹ = 0 Triangular superior
  9. Exemplos: −2 0 0 0 3 0 0 0 6

    Diagonal 3x3 3 0 0 0 0 3 0 0 1 8 4 2 2 0 5 1 Triangular inferior 4x4 1 4 2 5 2 0 3 1 0 7 0 0 0 0 0 0 5 0 0 9 2 0 3 1 4 Triangular superior 5x5
  10. Continuando no zoológico: Quando é igual à sua transposta: =

    . Simétrica Quando satisfaz: ≠ 0 ⟹ = > ,=1 0 Definida positiva Quando sua transposta é a inversa: = −1 ⟺ = = Ortogonal
  11. 1. + = + 2. ()= 3. = 4. (

    )= 5. = Observações: • De 1 e 2 decorre que matrizes simétricas constituem um subespaço do espaço de todas as matrizes. • A 5 significa que a matriz produto é uma matriz simétrica. Algumas propriedades da transposição:
  12. O produto de matrizes simétricas não é uma matriz simétrica!.

    E, além disso, também não é comutativo. Vejam os contra-exemplos! 3 2 2 1 0 1 1 −2 = 2 −1 1 0 0 1 1 −2 3 2 2 1 = 2 1 −1 0
  13. A matriz A = 1 0 0 0 −2 0

    0 0 3 não é definida-positiva, pois se = temos: = 1 0 0 0 −2 0 0 0 3 = = −2 3 = 2 − 22 + 32 e escolhendo = 1 2 0 obtemos = 0.
  14. A matriz = 5 −1 2 −1 2 2 2

    2 3 é definida-positiva, pois para = temos: = − 2 + + 2 2 + 2 + 2, uma soma de quadrados, que só se anula para = = 0 0 0 . Surfista, confira essa afirmação da Mestra!
  15. = sen cos 0 −cos sen 0 0 0 1

    sen − cos 0 cos sen 0 0 0 1 = 1 0 0 0 1 0 0 0 1 Loirinha, agora é tua vez de conferir a afirmação da Mestra! A matriz = sen − cos 0 cos sen 0 0 0 1 é ortogonal pois :
  16. Este programa recebe uma matriz A e um vetor b.

    Em seguida, chama solve(A,b) para resolver o sistema linear A*x = b.
  17. Mestra esse algoritmo para achar a inversa é válido para

    qualquer tipo de matriz, não? Sim minha filha e é o mais rápido! Mas você tem que usar o solve( ) para um sistema cuja matriz não possui características particulares.
  18. Uma matriz = [ ] possui banda inferior k quando

    = 0 para − > . E possui banda superior k, quando = 0 para − > . = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × Eis uma matriz A com banda inferior 2 (em laranja) e banda superior 1 (em verde):
  19. = × × 0 0 0 × × × 0

    0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ Se A é uma matriz 100x100 teríamos 10.000 elementos para armazenar. Nesse armazenamento de banda precisaremos apenas de 4x100-4 ∎ ∎ ∎ ∎ A ação para armazenar uma matriz de banda como essa é “esquecer os zeros e deixar cair”! Vejam só:
  20. Um zoológico bem complexo: A adjunta de = [ ],

    anotada AH, é definida por = [ ]. É a conjugada da transposta. Adjunta Quando satisfaz: = . É a correspondente de uma matriz simétrica, com elementos complexos, Hermitiana, ou auto-adjunta Quando sua adjunta é a inversa: = −1 ⟺ = = Unitária
  21. Álgebra linear Apresentação Rotinas básicas Decomposições Funções matriciais Matrizes especiais

    Agora vamos explorar as fatorações (Decompositions) existentes na scipy.linalg.
  22. Esta é uma das rotinas do scipy.linalg para realizar a

    fatoração PLU. Ela só realiza a fatoração; não apresenta a solução.
  23. A rotina pode efetuar o pivoteamento por linha (ou não).

    Ela retorna: • as três matrizes P, L e U • ou apenas duas, a PL que é a triangular inferior permutada e a U.
  24. Esta é uma rotina alternativa à anterior. Ela recebe a

    matriz A e retorna uma matriz contendo a fatoração LU e um vetor contendo os índices de pivoteamento.
  25. Esta rotina recebe a fatorada LU e o vetor piv

    fornecidos pela LU_factor( ). Recebe também o termo independente b e retorna a solução x.
  26. Carl Friedrich Gauss 30/04/1777 – 23/02/1855 Vamos encerrar este conjunto

    de transparências detalhando 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.
  27. 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. Primeiro só descreveremos o processo em linhas gerais, usando as equações.
  28. 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.
  29. 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 ~ .
  30. 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.
  31. 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, exceto na 1ª linha.
  32. O último passo consiste em resolver o sistema triangular. Algo

    extremamente simples. O processo consiste, portanto, em eliminarmos incógnitas nas equações. Daí o nome: processo de eliminação. Alguns textos falam em escalonamento.
  33. Sim! Como é feita, de fato, a eliminação? Mestres, vocês

    descreveram o processo em linhas gerais! Só papo...
  34. 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1

    ⋮ 2 ⋱ ⋮ ⋯ 1 2 ⋮ = 1 2 ⋮ A x = b Pois é... Agora vamos refazer tudo usando a forma matricial do sistema linear:
  35. 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 ⋱ ⋮ ⋮ ⋯
  36. | = 11 12 ⋯ 1 1 21 22 ⋯

    2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ Linha 2 O processo de eliminação consiste em efetuar as operações elementares sobre as linhas da matriz aumentada. Linha 1 Linha n
  37. As operações elementares são: 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
  38. 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.
  39. Quando pensamos nas equações, fica óbvio que as operações elementares

    não alteram a solução. Isto não é tão óbvio quando olhamos para a matriz aumentada.
  40. 1 /11 31 /11 Multiplicadores Observem que os multiplicadores foram

    escolhidos de forma a anular todos os 1 abaixo do 11 . O 11 é intitulado “pivot”. Essas operações exigem que o pivot seja não-nulo: 11 ≠ 0. Caso 11 = 0 trocamos a linha 1 com outra onde 1 ≠ 0.
  41. 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!
  42. 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:
  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. 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.
  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 ouvido falar disto!
  46. 1 /11 31 /11 Multiplicadores Assumam que o pivot 11

    = 10−3. Então ao dividirmos por 11 estaremos multiplicando por 1.000 (≅ 1024 = 210) os elementos 21 , 31 , … , 1 . 11 ≅ 2−10
  47. 1 /11 31 /11 Multiplicadores Em outras palavras, Mestre, os

    multiplicadores 21 /11 , 31 /11 , … , 1 /11 perderam 10 dos seus 23 preciosos bits de precisão. 11 ≅ 2−10
  48. Que horror Mestra. É por isso que os pivots tem

    que ser (homens) grandes! Sim minha Loirinha querida! Despreze pivots pequenos, eles estragam o jogo. Sempre vá em busca dos maiores possíveis.
  49. Uma forma de minimizar a propagação de erro consiste em

    localizar o 1 de maior valor absoluto e trocar a linha onde ele está com a linha 1. Esta é a estratégia conhecida na literatura como pivoteamento parcial por coluna. Na realidade estamos apenas apontando a “ponta do iceberg”. Analisaremos tudo isto em detalhes, mais adiante!
  50. Para facilitar o entendimento, vamos mostrar 1º a execução do

    programa. Depois discutiremos seu código. Estes são a matriz A e o vetor b de nosso exemplo.
  51. pivot Eis a 1ª etapa do processo de eliminação. As

    operações elementares visam colocar zeros na 1ª coluna, abaixo do pivot.
  52. Mestre, seu programa está errado! A matriz A atualizada deveria

    ter sido esta aqui, com zeros abaixo do pivot. Mil perdões Loirinha! Esqueci de dizer que usei o espaço dos zeros para guardar os multiplicadores.
  53. A 2ª etapa do processo de eliminação. As operações elementares

    visam colocar zeros na 2ª coluna, abaixo do pivot. Novamente, guardei os multiplicadores no lugar dos zeros.
  54. A etapa final do processo de eliminação. A solução x

    foi obtida por retro-substituição.
  55. Uma foto da parte final do processo: • A matriz

    A foi transformada na matriz triangular superior que eu marquei em azul. • E os multiplicadores estão em verde, abaixo da diagonal.
  56. Esta é a fase da eliminação. O loop em k

    define as etapas do processo de eliminação. É definido o pivot e, se ele é não-nulo, criamos uma matriz-coluna, Mult, para guardar os multiplicadores..
  57. No loop em j atualizamos cada elemento da linha i

    efetuando a operação elementar sobre ele. O loop em i percorre as linhas abaixo do pivot e define um multiplicador para cada linha.
  58. Ainda dentro do loop em i efetuamos a operação elementar

    sobre o elemento da linha i do termo independente b. E guardamos esse multiplicador na linha i da coluna do pivot da matriz A e também na linha i do vetor Mult.
  59. Agora a fase da retro-substituição. • Da última equação achamos

    diretamente o valor de x n . • Substituímos esse valor na penúltima equação e achamos o valor de x n-1 . • Assim por diante, de trás para frente, até chegar à 1ª equação, de onde achamos x 1 .