A SciPy.LinAlg

...

Transcript

1. A scipy.linalg Prof. Paulo R. G. Bordoni UFRJ
2. Aulas atrás, estudamos o básico da NumPy e da MatPlotLib. Iniciaremos agora nosso estudo da SciPy – Scientific Python.
3. E o que é, precisamente, Computação Científica Mestre? É o que você aprenderá ao longo deste curso, Loirinha. Vou começar pelo aspecto histórico!
4. http://history.siam.org/ A visão da Sociedade de Matemática Aplicada e Industrial – SIAM:
5. Procure por Scientific Computing na Web, Surfista. Farei a busca com o Google, Mestra.
6. Peguei o retorno do Google e separei o que achei interessante:
7. Comecei pela Wikipedia. Eis o que achei:
8. Nosso Mestre já citou Maple, Mathematica e MatLab. Como ele, a Wikipedia também indica Python e SciPy!
9. O 2º que separeifoi o site da INTEL de HPC. HPC é High Performance Computation
10. Separei algumas chamadas desse site que me interessaram – estão todas na mídia recente!
11. A SciPy é a 5ª da lista retornada pelo Google. E é reconhecida como uma biblioteca fundamental para computação científica! É com ela que trabalharemos Galileu. Atualíssima, veja a seguir:
12. A EuroSciPy 2016 foi na Alemanha e a de 2015 foi na Inglaterra.
13. A SciPy Latin America 2016 foi em Florianópolis de 16 à 20 de maio.
14. Este é o “site” da SciPy.
15. Clicando no ícone da SciPy Library, chegamos aqui. Agora escolha Documentation, Loirinha!
16. Surfista, em Documentation escolha o Guia de Referência da SciPy.
17. Todos os tópicos tratados no Tutorial, na última versão da SciPy:
18. Em Introduction temos uma breve descrição da SciPy. Lembrem-se, os Tutoriais são para novatos!
19. É importante assinalar as convenções de importação utilizadas por padrão.
20. E agora todos os tópicos tratados na Referência, na última versão da SciPy: Assinalei os sub- pacotes que estudaremos.
21. Antes de mergulharmos fundo nas matrizes, fale um pouco sobre Entrada/Saída Manuel. Vou listá-los, mas não entrarei em detalhes!
22. Outros métodos além dos padronizados para lidar com arquivos:
23. Métodos padronizados? Quais? Aqueles da Numpy referidos logo no início, filha:
24. Os métodos da Numpy para lidar com arquivos:
25. Continuação:
26. A parte final:
27. Passageiros com destino à scipy.linalg, embarque imediato pela plataforma 9 ¼.
28. Precisaremos pesquisar muito para descobrir tudo sobre Álgebra Linear Computacional na biblioteca do castelo.
29. Vou buscar com Google, Galileu! Surfista, antes disso, para chegarmos preparados, faça uma busca na Internet por Álgebra Linear Computacional.
30. Quando busquei por Numerical Linear Algebra, o 1º resultado foi da Wikipedia:
31. Não deixem de olhar as listas indicadas de software de análise numérica e de bibliotecas numéricas. Muita indicações importantes nesta dica!
32. As primeiras estações da ferrovia, antes de Hogwarts, são paradas obrigatórias, chamadas BLAS, LAPACK e NETLIB. Temos que desembarcar para conhecê-las.
33. Antes quero mostrar a Estação da Luz e a Júlio Prestes, na cidade de São Paulo. Hoje também são a sede do Museu da Língua Portuguesa e da Orquestra Sinfônica do Estado de São Paulo, OSESP.
35. A estação Leopoldina, no Rio. Se a Princesa visse o descaso iria chorar! Os responsáveis públicos pelo abandono deveriam ser presos!
36. A Estação da Central do Brasil, também no Rio. Dias depois do comício do Jango, veio o golpe de 64!
37. As bibliotecas fundamentais sobre Álgebra Linear Computacional são a BLAS e LAPACK.
38. • 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.
39. Fui buscar BLAS no Google! Assinalei o endereço da netlib!
40. Este é o principal repositório de software científico que a comunidade científica mundial dispõe:
41. Respondendo antes de você perguntar, Loirinha:
42. A BLAS na www.netlib.org
43. FAQ é comigo! Inclusive já marquei minha 1ª pergunta:
44. Eis sua resposta, Loirinha:
45. O endereço para chegar aqui: http://www.netlib.org/blas/ No final, os links para todas as rotinas da BLAS!
46. Sim, Surfista. Clicando num link, somos encaminhados para a rotina correspondente na LAPACK. O código está em FORTRAN.
47. O código, criado pelo poderosíssimo feiticeiro Jack Dongarra em 1978. A última atualização é de 2011.
48. A LAPACK na Wikipedia:
49. A continuação da LAPACK na Wikipedia:
50. O site oficial da LAPACK, claro: www.netlib.org.
51. Agora vamos à biblioteca do Castelo de Hogwarts. É lá que encontraremos o Tutorial da scipy.linalg.
52. Marquei o motivo para usar a linalg da SciPy. Mas já tem linalg na NumPy. Por que repetir?
53. Antes que você pergunte, Surfista, vá ler o resto desta discussão lá no Tutorial!
54. A estrutura do tutorial Linear Algebra é a seguinte: Álgebra linear Apresentação Rotinas básicas Decomposições Funções matriciais Matrizes especiais
55. Os tópicos em rotinas básicas são : 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
56. Vamos começar pelos determinantes, que já vimos na aula passada! E indo à referencia do scipy.linalg, nas rotinas básicas, temos:
57. A resposta à tua pergunta está ao pé da página Loirinha: por outro algoritmo. Mas vocês já provaram que o cálculo de determinantes é inviável, na prática. Então, como o pacote scipy faz isto?
58. E o que é fatoração LU? Mais adiante exploraremos tudo o que há na Linalg sobre fatorações. Em particular sobre a fatoração LU.
59. Respondendo à tua pergunta Loirinha, escrevemos como um produto de duas matrizes triangulares e : = ∙ A só tem 1’ na diagonal. Assim det = det ∙ = det ∙ det = det()
60. Neste programa entramos com uma matriz A, via teclado, e em seguida chamamos a função det(A).
61. Voltando às rotinas básicas da Scipy, vejamos o que há no tópico “Calculando normas”. 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
62. A SciPy possibilita calcular diversas normas. Tanto normas de vetores como de matrizes.
63. Por enquanto, vamos examinar apenas as normas vetoriais. Veremos as normas matriciais mais adiante. Para elas precisaremos de um pouco mais de teoria. É importantíssimo observar que todas aquelas cujo parâmetro ord é negativo não são normas.
64. Surfista, apresente alguns contra-exemplos que justifiquem a última afirmação do Mestre! Quando ord = inf temos a norma do máximo. Para ord = 1 temos a norma da soma e para ord = 2 temos norma euclidiana. Elas são anotadas ∞ , 1 e 2 respectivamente. As outras são anotadas , ( = ).
65. Já provamos que 2 = 1 2 + 2 2 + ⋯ + 2 define uma norma. A verificação que 1 = 1 + 2 + ⋯ + é uma norma é óbvia. Da mesma forma, a comprovar que ∞ = max{ 1 , 2 , ⋯ , } define uma norma também é fácil.
66. Para prová-la precisamos da desigualdade de Hölder: , ≤ ′ com ′ = −1 . A dificuldade para provar que = |1 | + |2 | + ⋯ + | | é uma norma está na desigualdade triangular: + ≤ + . Ela é conhecida na literatura como desigualdade de Minkowski.
67. Observem que a desigualdade de Hölder desempenha um papel análogo ao da desigualdade de Cauchy-Schwarz. Não vamos provar nenhuma das duas. Porém, se você estiver interessado, Surfista, já sabe as palavras-chave para buscar as demonstrações.
68. Na referência do SciPy encontramos os detalhes dos parâmetros para usar a função norma:
69. Por enquanto só veremos normas de vetores. Atenção para a referência!
70. Um exemplo de utilização da função norm( ) para vetores.
71. Mestre, parece que os valores numéricos das normas estão decrescendo ao valor da ∞ . É verdade minha filha. Prove que para ∀ ∈ ℝ temos lim →∞ = ∞
72. Para cada norma vetorial em ∙ : ℝ → ℝ podemos associar um subconjunto do próprio ℝ. Estou me referindo às bolas abertas de raio > 0 com centro num ponto ∈ ℝ ∙ , , = ∈ ℝ . . − <
73. Para a norma euclidiana em ℝ2, é o conjunto dos pontos em ℝ2 satisfazendo: ∙ 2, , = = (, ) ∈ ℝ2 . . − 2 < , sendo = , o centro E como 2 = 2 + 2, temos − 2 = ( − )2 + ( − )2.
74. Que é o conjunto dos pontos interiores à circunferência 2 + 2 = 1, de centro na origem e raio 1. Portanto a bola unitária na ∙ 2 , com centro na origem, corresponde ao conjunto definido pelos pontos , ∈ ℝ2 cujas coordenadas satisfazem a desigualdade 2 + 2 < 1
75. Fiz um programa que desenha a circunferência e bola unitária em ℝ2. Vejam o resultado:
76. O programa que desenha as bolas unitárias na ∙ :
77. Bolas na ∙ , para diversos
78. Surfista, não vá jogar bola com 1 ou ∞ . Você pode machucar seu pé! É, são bolas que escapam da nossa intuição euclidiana.
79. Note que as bolas “crescem’ da ∙ 1 para a ∙ ∞
80. Agora vamos explorar o que há na scipy.linalg para resolução de sistemas lineares.
81. Vamos ao tópico “Resolvendo sistemas lineares” de Rotinas básicas. 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
82. No Tutorial da SciPy, clicando em “Resolvendo sistemas lineares”, obtemos:
83. As rotinas básicas no Guia de Referência do scipy.linalg para resolução de sistemas lineares.
84. 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)
85. Recordando as definições: 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
86. 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
87. = 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:
88. 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:
89. 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 é, obviamente 1 = 2 = ⋯ = 100 = 1/
90. 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.
91. 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
92. Mestra, fiz até este programa: Veja a substituição
93. 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 .
94. Esta é rotina da SciPy para resolver sistemas triangulares.
95. O detalhamento dos seus parâmetros.
96. Eis meu programa:
97. Examinando as informações sobre o comando solve( ).
98. Este programa recebe uma matriz A e um vetor b. Em seguida, chama solve(A,b) para resolver o sistema linear Ax = b.
99. Outra forma de obter informações do solve( )
100. Uma forma de calcular a inversa de uma matriz A é resolvendo a equação matricial = usando solve( )
101. 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):
102. = × × 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ó:
103. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ ∎ ∎ ∎ ∎ Mestre, prefiro pensar em “deitar” as diagonais.
104. Esta é rotina específica da scipy.linalg para resolver sistemas lineares cuja matriz é de banda.
105. Novamente, os detalhes sobre os parâmetros
106. Um programinha usando a solve_banded( )
107. Um outro programa usando a solve_banded( )
108. Uma execução do programa.
109. Um exemplo maior! Nesse caso A possui 72 = 100-28 zeros.
110. 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:
111. 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:
112. 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
113. 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.
114. 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!
115. = 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 :
116. 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
117. As matrizes ortogonais (e as unitárias, no caso complexo) desempenham um papel fundamental em Álgebra Linear Computacional. Aguardem!
118. Esta rotina é apropriada para matrizes hermitianas (simétricas no caso real) e positivo-definidas.
119. Detalhamento dos parâmetros.
120. 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.
121. 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
122. Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
123. Uma execução do programa. Se a matriz A não for definida positiva, o programa acusa. Vejam a seguir.
124. O 4º elemento da diagonal principal é o responsável por tornar a matriz A não positiva definida. A linalg acusa o fato! ERRO
125. 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, = .
126. Tchau, até a próxima aula!