Álgebra Linear Computacional I

Por efetuar

Transcript

1. Álgebra linear computacional I Prof. Paulo R. G. Bordoni UFRJ
2. Aula retrasada iniciamos nossa viagem às terras da scipy.linalg. Há muita coisa a ser vista! Passageiros com destino à scipy.linalg.
3. Surfista, procure na Internet por mapas da Álgebra Linear Computacional no mundo da Computação Científica. Vou buscar no Google, Mestre!
4. Comecei com scientific computation, Mestre:
5. Python e SciPy!
6. http://history.siam.org/ Já com history of scientific computation, achei:
7. Quando busquei por Numerical Linear Algebra, o 1º resultado foi da Wikipedia:
8. • 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.
9. Fui buscar a BLAS no Google!
10. Um repositório de software importantíssimo!
11. O que é a Netlib:
12. A BLAS na www.netlib.org
13. É isso aí, Loirinha:
14. O que é a BLAS:
15. A LAPACK na Wikipedia:
16. A continuação de LAPACK na Wikipedia:
17. O site oficial da LAPACK, claro: www.netlib.org.
18. 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 !
19. 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!
20. No Tutorial da SciPy, clicando em “Resolvendo sistemas lineares”, obtemos:
21. 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
22. As rotinas básicas no Guia de Referência do scipy.linalg para resolução de sistemas lineares.
23. 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)
24. 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
25. 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
26. = 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:
27. 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:
28. 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/
29. 11 1 = 1 21 1 31 1 41 1 +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 Surfista, resolva o sisteminha 4x4 aos meus pés, cuja matriz é triangular inferior.
30. 11 1 = 1 21 1 31 1 41 1 +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 Mestra, assumindo ≠ 0 para = 1,2,3,4, temos: 1ª eq. ⇒ 1 = 1 /11 2ª eq. ⇒ 2 = (2 − 21 1 )/22 3ª eq. ⇒ 3 = 3 − (31 1 + 32 2 ) /33 4ª eq. ⇒ 4 = 4 − (41 1 + 42 2 + 43 3 ) /44
31. Mestra, fiz até este programa: Veja a substituição
32. 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 .
33. Esta é rotina da SciPy para resolver sistemas triangulares.
34. O detalhamento dos seus parâmetros.
35. Eis meu programa:
36. Examinando as informações sobre o comando solve( ).
37. Este programa recebe uma matriz A e um vetor b. Em seguida, chama solve(A,b) para resolver o sistema linear Ax = b.
38. Outra forma de obter informações do solve( )
39. Uma forma de calcular a inversa de uma matriz A é resolvendo a equação matricial = usando solve( )
40. 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):
41. = × × 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ó:
42. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ ∎ ∎ ∎ ∎ Mestre, prefiro pensar em “deitar” as diagonais.
43. Esta é rotina específica da scipy.linalg para resolver sistemas lineares cuja matriz é de banda.
44. Novamente, os detalhes sobre os parâmetros
45. Um programa usando a solve_banded( )
46. Uma execução do programa.
47. Um exemplo maior! Nesse caso A possui 72 = 100-28 zeros.
48. 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:
49. 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:
50. 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
51. 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.
52. 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!
53. = 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 :
54. 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
55. As matrizes ortogonais (e as unitárias, no caso complexo) desempenham um papel fundamental em Álgebra Linear Computacional. Aguardem!
56. Esta rotina é apropriada para matrizes hermitianas (simétricas no caso real) e positivo-definidas.
57. Detalhamento dos parâmetros.
58. 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.
59. 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. O nome do programa é: sist_banda_sim_def_pos.py
60. Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
61. Uma execução do programa. Se a matriz A não for definida positiva, o programa acusa. Vejam a seguir.
62. O 4º elemento da diagonal principal é o responsável por tornar a matriz A não positiva definida. A linalg acusa o fato! ERRO
63. 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, = .
64. Há muita coisa a ser vista em nossa viagem à scipy.linalg. A paisagem é magnífica!
65. Á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.
66. A famosa fatoração LU.
67. Rotinas existentes na Referência da scipy.linalg para a fatoração LU.
68. 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 = permite calcular det().
69. Mestre, que outras vantagens Maria Loirinha 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!
70. 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 já vimos, o maior trabalho reside em fatorar (i.é, efetuar o processo de eliminação).
71. Os detalhes para usar a linalg.lu( ):
72. O que linalg.lu( ) retorna:
73. A rotina pode efetuar o pivotamento 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.
74. Conferindo a fatoração LU:
75. Um 1º exemplo em que não necessidade de trocas de linhas,
76. 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.
77. Uma outra forma de chamar a fatoração LU da SciPy. Como no 1º exemplo, não houve a necessidade de pivotamento.
78. Chamando a 2ª forma da fatoração LU da SciPy. Agora com pivotamento.
79. Esta é lu_factor( ) que efetua a fatoração . Ela recebe a matriz A e retorna uma matriz contendo a fatoração LU e um vetor P contendo os índices de pivoteamento.
80. Os parâmetros retornados. Esta rotina usa a LAPACK.
81. 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.
82. Resolvendo um sistema linear Ax = b via fatoração LU com pivotamento.
83. Ela é boa para resolver sistemas Ax=b com muitos b’s, mas a mesma A.
84. Vejam, entrei com uma matriz 4x4 e 2 termos independentes. Poderiam ser muitos: 10, 50, etc.
85. Uma brevíssima explicação do que é a fatoração de Cholesky.
86. É 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!
87. 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?
88. 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.
89. As rotinas do SciPy para a fatoração de Cholesky.
90. Surfista e Loirinha, deixaremos com vocês a tarefa de construir alguns programas usando essas rotinas de Cholesky. É uma dica p’rá prova colega!
91. A explicação para passagem de parâmetros.
92. Esta rotina só efetua a fatoração. Deve ser usada com a próxima.
93. Esta rotina resolve o sistema após fatorado pela rotina cho_factor( ).
94. Esta é específica para matrizes de banda. Só fatora.
95. Esta resolve o sistema linear, após a fatoração pela cholesky_banded( ).
96. Tchau, até a próxima aula!