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

Álgebra linear computacional II

Álgebra linear computacional II

Iniciamos falando sobre a BLAS e a LAPACK e passamos às rotinas básicas da sipy.linalg para solução de sistemas lineares (solve, solve_triang., solve_banded, etc). Depois estudamos as fatorações LU e Cholesky.

Paulo Bordoni

October 10, 2013
Tweet

More Decks by Paulo Bordoni

Other Decks in Education

Transcript

  1. Aula passada iniciamos nossa viagem às terras da scipy.linalg. Há

    muita coisa a ser vista! Passageiros com destino à scipy.linalg.
  2. Surfista, procure na Internet por mapas da Álgebra Linear Computacional

    no mundo da Computação Científica. Vou buscar no Google, Mestre!
  3. • BLAS: Basic Linear Algebra Subprograms, • LAPACK: Linear Algebra

    Package. Repetindo: BLAS e LAPACK, bibliotecas altamente otimizadas que implementam os algoritmos básicos de Álgebra Linear Computacional. Não deixem de olhar as listas indicadas de software de análise numérica e de bibliotecas numéricas.
  4. 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, descobri o caminho das pedras. Muito mais que isto, percebi claramente a importância das bibliotecas BLAS e LAPACK !
  5. 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!
  6. Então vamos focar 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
  7. As características particulares da matriz do sistema linear definirão nosso

    roteiro para estudar os métodos de resolução do sistema. Tipo da matriz A em Ax = b: 1. Triangular 2. Genérica 3. De banda 4. Simétrica e definida positiva (hermitiana)
  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. = 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 det = 120 (= 1 × 3 × 5 × 2 × 4) = −3 0 0 0 −2 0 0 0 2 7 4 2 5 0 5 1 det = 0 (= −3 × 0 × 5 × 1) O determinante de uma matriz triangular é o produto dos termos da diagonal:
  11. Se A é uma matriz ordem n, então det =

    det (). Um sistema linear = possui uma única solução quando, e apenas, quando det () ≠ 0. Os dois teoremas abaixo, são fatos sobejamente conhecidos sobre sistemas lineares e determinantes. Entretanto determinantes conduzem a conclusões erradas sobre sistemas lineares. Vejam a seguir:
  12. Entretanto para = 0.01 temos det = 0.01100 = 10−200,

    um legítimo ZERO em qualquer computador! Então concluiríamos que o sistema não tem solução. 0 ⋮ 0 ⋯ 0 ⋮ ⋱ 0 0 0 0 1 2 ⋮ 100 = 1 1 ⋮ 1 A solução do sistema linear acima é, óbviamente 1 = 2 = ⋯ = 100 = 1/
  13. 11 1 = 1 21 1 31 1 41 1

    +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 É muito simples, Mestra. Se os ≠ 0, da 1ª eq. tiramos 1 = 1 /11 . Substituindo esse valor na 2ª eq. temos 2 = (2 − 21 1 )/22 . Depois 3 = 3 − (31 1 + 32 2 ) /33 . Finalmente, 4 = 4 − (41 1 + 42 2 + 43 3 ) /44 Só para saborear, vamos resolver o sisteminha 4x4 abaixo, cuja matriz é triangular inferior.
  14. Surfista, observe no código que o número de multiplicações cresce

    de 1 a n como uma PA de razão 1. Logo, o total de multiplicações é ( + 1)/2. Um algoritmo 2 .
  15. Este programa recebe uma matriz A e um vetor b.

    Em seguida, chama solve(A,b) para resolver o sistema linear Ax = b.
  16. Uma forma de calcular a inversa de uma matriz A

    é resolvendo a equação = usando solve( )
  17. 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):
  18. = × × 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ó:
  19. = × × 0 0 0 × × × 0

    0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ ∎ ∎ ∎ ∎ Mestre, prefiro pensar em “deitar” as diagonais.
  20. Quando é igual à sua transposta: = . Simétrica Quando

    satisfaz: ≠ 0 ⟹ = > ,=1 0 Definida positiva Quando sua transposta é a inversa: = −1 ⟺ = = Ortogonal Continuando a visita ao zoológico matricial:
  21. 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. • Atenção: a 5. não permite concluir que a matriz produto é uma matriz simétrica. Algumas propriedades da transposição:
  22. De fato, o produto de matrizes simétricas não é uma

    matriz simétrica!. E, além disso, o produto 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
  23. 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 = 2 1 0 obtemos = 0.
  24. 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!
  25. = 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 :
  26. 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
  27. As matrizes ortogonais (e as unitárias, no caso complexo) desempenham

    um papel fundamental em Álgebra Linear Computacional. Aguardem!
  28. O exemplo de armazenamento do Tutorial da scipy.linalg. Armazena-se apenas

    a diagonal e a banda superior (ou inferior). Em vez de 36 = 6 × 6 elementos, apenas 15 = 3 × 6 − 3.
  29. Esta 1ª parte do programa recebe uma matriz A (formato

    de banda), um termo independente b e devolve a solução x usando a rotina para matrizes simétricas, positivo-definidas e de banda da scipy.linalg.
  30. Esta 2ª parte apresenta a matriz A no formato padrão

    n x n e confere a solução, , verificando se, de fato, = .
  31. Uma execução do programa. Se a matriz A não for

    definida positiva, o programa acusa. Vejam a seguir.
  32. O 4º elemento da diagonal principal é o responsável por

    tornar a matriz A não positiva definida. A linalg acusa o fato! ERRO
  33. Que péssimo exemplo Mestre! Por que não incluiu um tratamento

    de erro, para evitar que seu programa pipocasse assim? Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
  34. Há muita coisa a ser vista em nossa viagem à

    scipy.linalg. A paisagem é magnífica!
  35. Á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.
  36. Mestre, que vantagem Maria leva com essa tal de fatoração

    LU? Antes do Mestre te responder, jovem, proponho que discuta com seus colegas sobre essa cultura de levar vantagem em tudo. Grande Filósofo!
  37. Respondendo tua pergunta, Loirinha: Tendo a fatoração = , resolver

    o sistema linear = é equivalente a resolver = . Então resolvemos primeiro o sistema = e depois o = . Note que esses dois são sistemas triangulares, facílimos de resolver! Conforme veremos depois, o maior trabalho reside em achar a fatoração.
  38. Sim, de = , segue que det () = det

    det det (). Veremos que L só possui 1’s na diagonal e que P é ortogonal. Assim, det = ± det = ±11 22 ⋯ A fatoração = também permite calcular det ().
  39. 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.
  40. 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.
  41. Já neste exemplo é evidente a necessidade de permutar a

    L0 com a L2 e depois a L2 com a L3. Atente para as matrizes A e P.
  42. Uma outra forma de chamar a fatoração LU da SciPy.

    Como no 1º exemplo, não houve a necessidade de pivotamento.
  43. 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.
  44. 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.
  45. É essa a ideia Surfista. Além disso, vale a volta;

    veja na próxima transparência! Confirmando, se A é simétrica e positiva definida está garantida a existência de uma matriz triangular superior tal que = . É como se U fosse a raiz quadrada de A, Mestre!
  46. Suponha que uma matriz não-singular A possa ser fatorada na

    forma = , com U triangular superior. Então A será simétrica e definida positiva. É difícil provar isto Mestra?
  47. E também que A é definida positiva, isto é ≠

    0 ⟹ > 0. De fato assumindo ≠ 0 e = teremos ≠ 0, devido à não-singularidade de U. Então = () = = = 2 > 0. Bem, minha querida, de = segue que: • A é simétrica, pois = ()= ()= = • U é não-singular, pois det ()2 = det() = det () ≠ 0.
  48. Surfista e Loirinha, deixaremos com vocês a tarefa de construir

    alguns programas usando essas rotinas de Cholesky. É uma dica p’rá prova colega!