Álgebra linear computacional I

Iniciamos mostrando a impossibilidade de resolver sistemas lineares pela regra de Crammer. Depois apresentamos a biblioteca SciPy e uma introdução à sicpy.linalg. Apresentamos também a BLAS e a Linpack.

Transcript

1. Álgebra linear computacional Prof. Paulo R. G. Bordoni UFRJ
2. Na aula passada começamos a olhar para matrizes. Vimos também como criar matrizes e efetuar operações matriciais básicas com NumPy. Vamos adiante!
3. Cramer apresentou sua fórmula para resolver sistemas lineares em 1750. Antes disso, Cardano, no final do século XVI, trabalhou com determinantes 2x2. Leibniz usou determinantes em 1693 para resolver sistemas lineares.
4. Um determinante de 2ª ordem, é denotado e definido por: Notem que o cálculo exige 2 multiplicações e 1 subtração. No computador serão 3 flops.
5. ≠ A regra de Cramer para resolver um sistema 2x2 como este é (prove, Surfista): O cálculo exige 3x2 multiplicações, 3 adições e 2 divisões. Ao todo 11 flops.
6. = 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 + = + = − obtemos:
7. Um determinante de 3ª ordem, é definido por: Desenvolvendo os determinantes de 2ª ordem, obtemos:
8. A regra de Cramer para resolver um sistema 3x3 como o acima é dada por:
9. Portanto para resolver um sistema 3x3 pela regra de Cramer são necessárias: (3/2)4! multiplicações, 4!-4 adições e 3 divisões. Notem o fatorial no número de flops. Foram necessárias (3/2)3! multiplicações e 3!-1 adições/subtrações para cada determinante 3x3.
10. Este é um sistema linear constituído por n equações lineares à n incógnitas. Sendo breve: um sistema n x n
11. Ax = b Aliás, Lewis Carrol conta num livro, que a notação para matrizes que usamos hoje foi proposta por Leibniz. Na forma matricial esse sistema se reescreve:
12. O determinante det A de uma matriz A de ordem n é definido recursivamente, como na próxima transparência.
13. Para n > 1, o desenvolvimento pela − é linha do det é dado por det = 1 1 + 2 2 + ⋯ + onde = (−1)+ e é o determinante de ordem n-1 obtido a partir da matriz A excluindo-se sua − é linha e − é coluna. Há um desenvolvimento semelhante a partir da − é coluna. Para n = 1, det A = a 11 .
14. det (ℴ 4 ) α det (ℴ 3 ) β det (ℴ 3 ) γ det (ℴ 3 ) δ det (ℴ 3 ) + − + − + − α det (ℴ 2 ) β det (ℴ 2 ) γ det (ℴ 2 ) det (ℴ 3 ) + − + − α det (ℴ 1 ) β det (ℴ 1 ) det (ℴ 2 ) + − 2 mult. 1 adiç. 3 mult. 2 adiç. 4 mult. 3 adiç. O cálculo de um determinante de ordem 4 envolve: 4! multiplicações 3! adições/subtrações
15. A regra de Cramer para resolver um sistema linear n x n como este é: Se = det ≠ 0 a solução é dada por 1 = 1 , 2 = 2 , ⋯ , = onde é o determinante da matriz A com a coluna k subtituída pelo termo independente b.
16. Portanto, para resolver um sistema linear n x n precisaremos efetuar o cálculo de n+1 determinantes de ordem n e n divisões. O total de operações envolvidas é então: + 1 ! multiplicações + 1 ( − 1)! adições n divisões
17. Bobagem! Para o Titan é mole! O tempo gasto, por qualquer computador, para resolver sistemas 30 x 30 é inimaginável !
18. A tabela a abaixo mostra o tempo gasto por um computador, que faz uma multiplicação em 10-9 segundos, para resolver um sistema n x n: Para um sistema 25x25 o tempo é de 500 milhões de anos. Comparem com a idade do universo! n tempo 10 0.004 seg. 15 22 min. 20 77 anos 25 0,5x109 anos
19. Definitivamente, é inviável resolver sistemas lineares pela regra de Crammer! Jamais esqueça de contar o número de operações envolvidas num algoritmo!
20. Nesta e nas próximas aulas mergulharemos de cabeça na Álgebra linear computacional. Utilizaremos a SciPy.
21. SciPy.org é o site da SciPy que é um codinome para Scientific Python.
22. SciPy é uma biblioteca construída a partir de NumPy. Já vem instalada com Python(x,y).
23. Na página de abertura da SciPy.org uma homenagem merecida ao criador da MatPlotLib.
24. Palestras no congresso, agora em 2013.
25. Mais palestras do congresso.
26. Um dos tutoriais do Congresso.
27. O Guia de referência da SciPy.
28. Todos os tópicos tratados na última versão da SciPy:
29. Clicando em Tutorial somos conduzidos a uma repetição do seu índice. De imediato vamos olhar para Introduction
30. Uma breve descrição de SciPy. Convenções de importação.
31. Em nosso curso trabalharemos com os sub pacotes que marquei.
32. É hora de embarcar no pacote scipy.linalg. Ele possui a álgebra linear computacional tão anunciada.
33. Tudo começa aqui. Já, já entrarei em mais detalhes sobre as BLAS e LAPACK.
34. Eis o motivo para usar a linalg da SciPy. Mas já tem linalg na NumPy. Por quê repetir?
35. Antes que você pergunte, Surfista, vá ler o resto desta discussão lá no Tutorial!
36. A estrutura do tutorial Linear Algebra, (scipy.linalg), é a seguinte: Álgebra linear Apresentação Rotinas básicas Decomposições Funções matriciais Matrizes especiais
37. 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 Vamos começar pelos determinantes.
38. O texto é um resumo do que já vimos no início da aula.
39. Então, indo ao Reference Guide do pacote Linear Algebra, em Basics, vamos encontrar os detalhes para usar a função. Mas se o cálculo de determinantes é inviável, na prática, como o pacote scipy faz isto?
40. A utilização da função det( ) é trivial e a resposta à sua pergunta está marcada no pé da página, Loirinha.
41. E o que é fatoração LU? É algo tão importante e fundamental em Álgebra linear computacional como respirar é para a vida! O Mestre informou que será o assunto da próxima aula.
42. Neste programa entramos com uma matriz A, via teclado, e em seguida chamamos a função det(A).
43. Voltando às rotinas básicas do 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
44. O SciPy possibilita calcular diversas normas. Tanto normas de vetores como de matrizes.
45. 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.
46. Surfista, apresente alguns contra-exemplos que justifiquem ú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 , ( = ).
47. 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.
48. 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.
49. 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.
50. Na referência do SciPy encontramos os detalhes dos parâmetros para usar a função norma:
51. Por enquanto só veremos normas de vetores. Atenção para a referência!
52. Um exemplo de utilização da função norm( ) para vetores.
53. Mestre, parece que os valores numéricos das normas estão tendendo ao valor da ∞ ! É verdade minha filha. Prove que para ∀ ∈ ℝ temos lim →∞ = ∞
54. 1 y x A bola unitária no espaço euclidiano ℝ é o conjunto 2 = ∈ ℝ | 2 ≤ 1 . Ela descreve o conjunto 2 + 2 + 2 ≤ 1, que nossa intuição entende por uma bola (de futebol) no espaço euclidiano ℝ3. No plano seria como um CD (dos Beatles).
55. 1 1 -1 -1 y x A bola unitária no espaço ℝ com a norma da soma é o conjunto 1 = ∈ ℝ | 1 ≤ 1 . No ℝ2 ela corresponde ao conjunto definido pela desigualdade + ≤ 1, mostrado na figura.
56. 1 1 -1 -1 y x A bola unitária no espaço ℝ com a norma do máximo é o conjunto ∞ = ∈ ℝ | ∞ ≤ 1 . No ℝ2 ela corresponde ao conjunto definido pela desigualdade max { , } ≤ 1, mostrado na figura.
57. Surfista, não vá jogar bola com 1 ou ∞ . Você pode se machucar! É, são bolas que escapam da nossa intuição euclidiana.
58. Este livro é o clássico dos clássicos sobre álgebra linear computacional. Vamos até reproduzir seu Sumário.
60. Este é outro clássico peso-pesado em álgebra linear computacional. Também vamos reproduzir seu Sumário.
61. São livros de cabeceira para quem for fazer pós-graduação em computação científica, matemática aplicada ou análise numérica.
62. O conteúdo desses dois livros é a definição mais apropriada de Álgebra linear computacional. Invista em você. Adquira os dois livros!
63. Um pouco mais de teoria. Vamos recordar conceito de transformação linear.
64. O conceito de função é, com certeza, o mais importante da matemática. Elas são o objeto do Cálculo diferencial e integral. Em particular, vamos considerar funções de um espaço vetorial U em um outro espaço vetorial V : : → , ∈ ⟼ = () ∈ .
65. Em outras palavras, funções : → tais que: • + ⟼ ( + ) = + , • ⟼ = (). Tais funções recebem o nome especial de transformações lineares. Particularizando ainda mais: funções : → de um espaço vetorial U para outro V que preservam as operações nativas de U e V.
66. = A ℝ ℝ Matrizes × definem transformações lineares de ℝ em ℝ através da multiplicação matriz x vetor. Uma matriz 3 x 4 define uma transformação linear ∶ ℝ4 → ℝ3, através da multiplicação: 1 2 3 = 11 12 13 14 21 22 23 24 31 32 33 34 1 2 3 4
67. Lembrem-se que: Se A é uma matriz m x n então: • + = + , ∀, ∈ ℝ • = , ∀ ∈ ℝ , ∀ ∈ ℝ Que são as condições de linearidade.
68. Um exemplo simples de ℝ2 em ℝ2 é 1 2 = −1 2 2 −3 1 2 que corresponde ao par de igualdades 1 = −1 + 22 2 = 21 − 32
69. − y x Em ℝ2 e em ℝ3 podemos visualizar o efeito de uma matriz como transformação linear. Por exemplo, a matriz = 1 0 0 −1 faz uma reflexão no eixo-x, pois − = 1 0 0 −1
70. O programa da próxima transparência permitirá a visualização da ação de uma matriz M sobre um triângulo ABC. Mostraremos alguns exemplos com diferentes matrizes M. Surfista, faça outros exemplos exploratórios.
71. O programa a seguir mostra a ação de uma matriz M de ordem 2 sobre um triângulo ABC, retângulo em B.
72. Em todos os exemplos estaremos usando a mesma escala nos dois eixos. A matriz M escolhida faz uma reflexão no eixo-y.
73. A matriz M escolhida faz uma contração no eixo-x e uma dilatação no eixo-y.
74. Já a ação desta matriz M é complicada de descrever.
75. Surfista, ajude a Loirinha. Vá ao programa que mostra a ação da matriz M sobre o triângulo ABC e marque os vértices correspondentes com o mesmo símbolo. Mestre, seu programa não mostra a correspondência entre os vértices originais e suas imagens!
76. Nas próximas três transparências comentaremos alguns detalhes do programa feito pelo Mestre. São detalhes técnicos...
77. A função plot( ) só funciona com arrays; não produz gráficos com matrizes. Vejam:
78. Face ao exposto na transparência anterior, precisamos transformar M_Triang[0] e M_Triang[1] em arrays, usando ravel( ).
79. A função fill( ) da MatPlotLib preenche uma poligonal fechada com a cor definida em rgb pela “string” hexadecimal ‘#rrggbb’. Assinalei um truque “sujo” do Mestre para centralizar o gráfico. Troque ‘.w’ por ‘+r’ e execute novamente para entender!
80. x x I A matriz identidade I é uma matriz quadrada de ordem n com 1’s na diagonal e 0’s fora dela. A identidade, como transformação linear não faz nada! ℝ ℝ I x = x
81. Uma matriz M muito próxima da identidade = 1 0 0 1 . Se M fosse a identidade I veríamos apenas o triângulo rosa, cobrindo o azul.
82. x y A z B ℝ ℝ ℝ BA O produto BA de uma matriz B por uma matriz A corresponde à aplicação linear composta de A seguida por B:
83. Lembrem-se que a multiplicação de matrizes não é comutativa. Na aula passada exibimos um programa que mostra que, usualmente, ≠ . Por exemplo, para = 1 2 1 2 e = −1 1 1 −1 temos = 1 −1 1 −1 e = 0 0 0 0 .
84. A inversa de uma matriz quadrada A de ordem n é uma matriz quadrada B de ordem n tal que ∙ = ∙ = Existem matrizes que não possuem inversa – são não-inversíveis. Uma matriz não possui mais que uma inversa. Além disso, a inversa da inversa é a própria!
85. Prova-se que A é inversível quando, e apenas quando, é não-singular, i.é, det () ≠ 0 Entretanto, verificar se det () ≠ 0 para saber se A é inversível é uma técnica nunca utilizada em Álgebra linear computacional.
86. x y A A-1 ℝ ℝ x y A x A-1 ℝ ℝ ℝ I Pela própria definição, a inversa A-1 de uma matriz A , quando considerada como aplicação linear, desfaz a ação de A.
87. Voltando às rotinas básicas do foquemos nossa atenção no tópico “Achando a inversa”. 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
88. Vejam o que Tutorial do scipy.linalg informa sobre a inversa de uma matriz:
89. Indo à Referência e clicando em inv( ) obtive as informações abaixo:
90. Obtendo a inversa de uma matriz A e conferindo.
91. Pedi ao Mestre para fazer um programa parecido àquele que mostra o efeito visual da ação de uma matriz num triângulo. Mas para a matriz inversa! Ele fará Loirinha, você é a queridinha dele!
92. O seu pedido satisfeito, Loirinha!
93. Verificando as ações de uma matriz e sua inversa, como transformações lineares. O deslocamento é um artifício para enxergarmos o efeito da inversa (em verde).
94. Surfista, descubra pela Internet o caminho das pedras para Computação Científica e Álgebra Linear Computacional. Vou buscar no Google, Mestre!
95. Com a palavra-chave scientific computation, achei:
96. Python e SciPy!
97. http://history.siam.org/ E com history of scientific computation, achei:
98. Quando busquei por Numerical Linear Algebra, o 1º resultado foi da Wikipedia:
99. • 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.
100. Fui buscar a BLAS no Google!
101. Um repositório de software importantíssimo!
102. O que é a Netlib:
103. A BLAS na www.netlib.org
104. Perguntas básicas:
105. O que é a BLAS:
106. A LAPACK na Wikipedia:
107. A continuação de LAPACK na Wikipedia:
108. O site oficial do LAPACK, claro: www.netlib.org
112. Tchau, até a próxima aula!