Introdução à SciPy Linalg

...

Transcript

1. Introdução à SciPy Prof. Paulo R. G. Bordoni UFRJ
2. Antes da SciPy vou mostrar como usamos a NumPy carregar arquivos do disco para a RAM e como salvar da RAM para o disco.
3. Os métodos disponibilizados pela NumPy para carregar e salvar arquivos:
4. Vou começar pelo load( ): Mostre-nos como usá-los, Mestre!
5. Esta é a rotina load( ):
6. A continuação:
7. E aqui está a rotina save( ):
8. Mostre-nos um código que usa as duas rotinas, Mestre Faremos um programa que: • recebe via teclado um vetor e uma matriz e, pequenos; • Salva no disco rígido e em arquivos distintos de nome arq_vetor e arq_matriz ; • Carrega e do disco rígido para a RAM, mas com os nomes e ; • Mostra e para comparar com e .
9. O código que você pediu, Surfista. O código ficou longo porquê ela caprichou na didática!
10. Veja, Surfista, uma execução do programa. Agora, repita o programa, sem os aspectos didáticos!
11. Boa pergunta, Loirinha! Nesse caso você usaria uma das rotinas • savez( ) • savez_compressed( ). Confira no código a seguir: Mestre, e seu quisesse guardar o vetor e a matriz num arquivo só?
12. O código que você pediu, Loirinha. Novamente, o código ficou longo porquê ela caprichou na didática!
13. Veja, Loirinha, uma execução do programa. Agora, repita o programa, sem os aspectos didáticos!
14. Vamos mergulhar na SciPy, Scientific Python:
15. 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!
16. http://history.siam.org/ A resposta da Sociedade de Matemática Aplicada e Industrial – SIAM:
17. Outra resposta, focando em programação:
18. A continuação:
19. Conferências da SciPy pelo mundo:
20. Este é o “site” da SciPy.
21. Clicando no ícone da SciPy Library, chegamos aqui. Agora escolha Documentation, Loirinha!
22. Surfista, em Documentation escolha o Guia de Referência da SciPy.
23. Todos os tópicos tratados no Tutorial, na última versão da SciPy:
24. Atenção para as as convenções de importação:
25. Na referência marquei em vermelho os sub-pacotes que estudaremos.
26. Passageiros com destino à scipy.linalg, embarque imediato pela plataforma 9 3/4 .
27. Precisaremos pesquisar muito para descobrir tudo sobre Álgebra Linear Computacional na biblioteca do castelo.
28. Vou buscar com Google, Galileu! Surfista, antes disso, para chegarmos preparados, faça uma busca na Internet por Álgebra Linear Computacional.
29. Quando busquei por Numerical Linear Algebra, o 1º resultado foi da Wikipedia:
30. 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!
31. As primeiras estações da ferrovia, antes de Hogwarts, são paradas obrigatórias, chamadas BLAS e LAPACK . Foram marcadas em verde. Temos que desembarcar para vê-las.
32. 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.
34. A estação Leopoldina, no Rio. Se a Princesa visse o descaso iria chorar! Os responsáveis públicos pelo abandono deveriam ser presos!
35. A Estação da Central do Brasil, também no Rio. Dias depois do comício do Jango, veio o golpe de 64!
36. As bibliotecas fundamentais sobre Álgebra Linear Computacional são a BLAS e LAPACK.
37. • BLAS: Basic Linear Algebra Subprograms, • LAPACK: Linear Algebra Package, • www.netlib.org Repetindo: BLAS e LAPACK, bibliotecas altamente otimizadas que implementam os algoritmos básicos de Álgebra Linear Computacional. Estão na netlib.
38. Este é o principal repositório de software científico que a comunidade científica mundial dispõe:
39. Respondendo antes de você perguntar, Loirinha:
40. A BLAS na www.netlib.org
41. FAQ é comigo! Inclusive já marquei minha 1ª pergunta:
42. Eis sua resposta, Loirinha:
43. A LAPACK na Wikipedia:
44. A continuação da LAPACK na Wikipedia:
45. Agora vamos à biblioteca do Castelo de Hogwarts. Vamos descobrir tudo sobre a scipy.linalg.
46. Passei antes pela Biblioteca Nacional, na Av. Rio Branco, no centro do Rio de Janeiro.
47. Ela é grandiosa. Monumental !!
48. Por esse motivo, ela é enorme. São 92 rotinas ao todo, mas vou listá-las ! A linalg é a biblioteca da scipy sobre Álgebra Linear.
49. As rotinas básicas no Guia de Referência da scipy.linalg.
50. Continuação
51. Problemas de auto- valores e auto-vetores
52. Fatorações ou decomposições
53. Continuação
54. Funções matriciais
55. Resolução de equações matriciais
56. Algumas matrizes especiais
57. Onde buscar as rotinas da blas e da lapack na scipy.linalg .
58. Procure por scipy.linalg.blas
59. Procure por scipy.linalg.lapack
60. Vamos começar pelas rotinas mais simples.
61. Repetindo: pelas mais simples. Sim. Pelo cálculo de normas de vetores.
62. = 1 2 1 2 Para um vetor = 1 2 ela é anotada 2 calculada através do teorema de Pitágoras 2 = 1 2 + 2 2 A norma euclidiana é a forma natural de calcular o tamanho de um vetor.
63. A generalização para um vetor = 1 2 … ∈ ℝ é óbvia: 2 = 1 2 + 2 2 + ⋯ + 2 Outras normas importantes são: • a norma da soma, anotada 1 • a norma do máximo, anotada ∞
64. Para um vetor = 1 2 … ∈ ℝ essas normas são calculadas de forma bem mais simples que a 2 : 1 = 1 + 2 + ⋯ + Norma da soma ∞ = 1 , 2 , , ⋯ , Norma do máximo
65. Além dessas três existem as , para ∈ ℕ, que generalizam a norma euclideana, e são calculadas como = 1 + 2 + ⋯ + Existem também as importantíssimas normas de energia, que utilizam formas quadráticas definidas através de matrizes.
66. Para provar que todas elas são normas é preciso garantir as quatro condições abaixo: Uma função : → ℝ, definida num espaço vetorial é uma norma quando satisfaz as condições: a) ≥ 0, ∈ b) = 0 ⇒ = 0 , ∈ c) ∙ = ∙ , α ∈ ℝ e ∈ d) + ≤ + , , ∈
67. As funções norma na scipy.linalg.
68. Por enquanto veremos apenas as normas vetoriais.
69. Um exemplo de utilização da função norm( ) para vetores.
70. Mestre, parece que os valores numéricos das normas estão decrescendo ao valor da ∞ . É verdade minha filha. Prove que para ∀ ∈ ℝ temos lim →∞ = ∞
71. 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 ∈ ℝ ∙ , , = ∈ ℝ . . − <
72. 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.
73. 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
74. Fiz um programa que desenha a circunferência e bola unitária em ℝ2. Vejam o resultado:
75. O programa que desenha as bolas unitárias nas ∙ :
76. Bolas na ∙ , para diversos
77. Surfista, não vá jogar bola com 1 ou ∞ . Você pode machucar seu pé! É, são bolas que escapam da nossa intuição euclidiana.
78. Note que as bolas “crescem’ da ∙ 1 para a ∙ ∞
79. Como já disse, as bolas “crescem’ da ∙ 1 para a ∙ ∞ mas o tamanho de um vetor decresce da ∙ 1 para a ∙ ∞ . Convença-se, Surfista, que esses resultados não são contraditórios!
80. Agora, vamos começar a estudar algoritmos para a resolução de sistemas de equações lineares.
81. A todo sistema linear está associada uma matriz, a matriz dos coeficientes do sistema, e dois vetores coluna, o termo independente e vetor das incógnitas. 11 1 + 12 2 + ⋯ 1 = 1 21 1 + 22 2 + ⋯ 2 = 2 31 1 + 32 2 + ⋯ 3 ⋯ 1 1 + 2 2 + ⋯ = ⋮ = 3 ⋯
82. 11 1 + 12 2 + ⋯ 1 = 1 21 1 + 22 2 + ⋯ 2 = 2 31 1 + 32 2 + ⋯ 3 ⋯ 1 1 + 2 2 + ⋯ = ⋮ = 3 ⋯ No caso do sistema acima, temos: = 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 ⋮ 2 ⋱ ⋯ ⋮ Matriz dos coeficientes = 1 2 ⋮ Termo independente = 1 2 ⋮ Vetor das inógnitas
83. Vamos começar pela regra de Cramer. A regra de Cramer é um algoritmo para resolver sistemas de lineares que utiliza determinantes.
84. Cramer apresentou sua fórmula para resolver sistemas lineares em 1750. Hoje conhecida como “a regra de Cramer” Antes disso, Cardano, no final do século XVI, trabalhou com determinantes 2x2. Leibniz usou determinantes em 1693 para resolver sistemas lineares.
85. Portanto, para utilizar a regra de Cramer precisamos saber como calcular determinantes. No caso de matrizes 2x2 é imediato: Se = 11 12 21 22 então det = 11 22 − 21 12 .
86. Se = 11 12 13 21 22 23 31 32 33 , uma matriz 3x3, usamos recursividade, isto é, recaímos no cálculo de determinantes 2x2. Mostre o “desenvolvimento pela 1ª linha”, Mestre. Ei-lo: det = 11 det 11 − 12 det 12 + 13 det 13 onde 11 = 22 23 32 33 , 12 = 21 23 31 33 e 13 = 21 22 31 32
87. Juro que não entendi como esses 11, 12, 13 aparecem, Mestre! Fácil, 1 é a matriz que sobra da matriz quando excluímos sua linha 1 e coluna j . Confira com 12: = 11 12 13 21 22 23 31 32 33 ⇒ 21 23 31 33 = 12
88. E o sinal de −, como apareceu ? O sinal de cada parcela é dado por (−1)1+
89. Então, para uma matriz 4x4 sigo a mesma ideia ? Garoto esperto, Exatamente!
90. Conferindo, se = 11 12 13 21 22 23 31 41 32 42 33 43 41 42 43 44 então: det = 11 det 11 − 12 det 12 + 13 det 13 − 14 det(14) Só que agora 11, 12, 13, 14 são matrizes 3x3.
91. Ah, Galileu é o computador que fará as contas. É rapidinho!!! Confira a quantidade de contas que você fará, Surfista.
92. É: 1 = 1 , 2 = 2 para ≠ 0. A regra de Cramer para o sistema 2x2 ቊ 11 1 + 12 2 = 1 21 1 + 22 2 = 2 Onde: = det 11 12 21 22 , 1 = det 1 12 2 22 , 2 = det 11 1 21 2
93. = 2 3 1 −1 2 1 1 3 = −5 5 = −1 = 3 1 −1 3 2 1 1 3 = 10 5 = 2 Aplicando a regra de Cramer para resolver o sistema 2x2 ቊ 2 + = 3 + 3 = −1 obtemos:
94. 1 = 1 , 2 = 2 , 3 = 3 , ≠ 0 A regra de Cramer para resolver o sistema ቐ 11 1 + 12 2 + 12 3 = 1 21 1 + 22 2 + 23 3 = 2 31 1 + 32 2 + 33 3 = 3 é: Só que agora = det onde é a matriz dos coeficientes do sistema linear: = 11 12 13 21 22 23 31 32 33
95. 1 = 1 12 13 2 22 23 3 32 33 , 2 = 11 1 13 21 2 23 31 3 33 , 3 = 11 12 21 22 2 31 32 3 E quem são 1 , 2 , 3 ? São os determinantes das matrizes 1 , 2 , 3 obtidas substituindo as colunas 1, 2, 3 pelo termo independente do sistema:
96. No braço dará muito trabalho. Mas, feito o código, no computador é moleza, tudo ocorrerá num milionésimo de piscar de olhos! Insisto que você calcule a quantidade de “continhas”, Surfista.
97. 1 = 1 , 2 = 2 , 3 = 3 , 4 = 4 , ≠ 0 Para resolver o sistema 4x4, 11 1 + 12 2 + 12 3 + 14 4 = 1 21 1 + 22 2 + 23 3 + 24 4 = 2 31 1 + 32 2 + 33 3 + 34 4 = 3 41 1 + 42 2 + 43 3 + 44 4 = 4 por Cramer, tudo se repete: Só que agora com matrizes 4x4.
98. Vou mostrar para vocês porque o Galileu está insistindo tanto para que vocês contem a quantidade de operações (as continhas...).
99. Sim. Confiram que para um determinante 2 × 2 fazemos 2 multiplicações. E para o determinante de uma matriz 3 × 3, pela fórmula det = 11 det 11 − 12 det 12 + 13 det 13 , fazemos 6 = 3 × 2 multiplicações.
100. Pensando na fórmula de recorrência para calcular o determinante de uma matriz 4 × 4 é imediato que faremos 24 = 4 × 6 multiplicações, pois serão 4 determinantes 3 × 3. = 11 11 − 12 12 + 13 13 − 14 det(14)
101. Portanto, você fará 24 = 4 × 6 = 4 × 3 × 2 = 4! multiplicações. Para o determinante de uma matriz de ordem , serão ! multiplicações. Claro Sherlock! Mas as contas serão realizadas pelo computador.
102. Fiz este programinha para apoiar a argumentação do Galileu!
103. O gráfico em vermelho, aos meus pés, é o da função Γ: ℝ+ → ℝ cuja restrição ao conjunto ℕ dos naturais é a função fatorial.
104. Uma tabela dessas funções para inteiro, variando de 0 à 30. Só para contar o número de dígitos
105. Surfista, contei a quantidade de dígitos apenas das multiplicações para calcular o determinante de matrizes 20 × 20, 25 × 25, 30 × 30: Veja, são números enormes!!! 20 × 20 → 20! → 19 í 25 × 25 → 25! → 26 í 30 × 30 → 30! → 33 í
106. Mas o Titan – Cray XK7 opera em teraflop/s vai tirar isto de letra!
107. Este programa fornece o tempo para o Titan calcular os determinantes:
108. Repetindo: O tempo gasto pelo Cray XK7, para calcular o determinante de uma matriz n x n na expansão por linhas é: n tempo 20 48,41 min. 25 7,37 séculos 30 1,076 idade univ.
109. Jamais esqueça de contar o número de operações efetuadas num algoritmo!
110. Vamos examinar como a SciPy permite calcular determinantes. O Galileu acabou de provar a inviabilidade do cálculo ! Como a scipy calcula determinantes, Mestra ?
111. Por outro algoritmo, Loirinha, indiquei ao pé da página.
112. E o que é fatoração LU? Assunto das próximas aulas, Loirinha.
113. Adiantando informação, Loirinha, fatoramos como um produto de duas matrizes triangulares e : = ∙ . A matriz L é triangular inferior (o L é de lower) e a U triangular superior (o U é de upper). A só tem 1’ na diagonal. Assim det = det ∙ = = det det = = det = 11 22 ⋯ .
114. Neste programa entramos com uma matriz A, via teclado, e em seguida chamamos a função det(A).
115. Tchau, até a próxima aula!